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SUMMARY 


The  IAR  Twin  Otter  Atmospheric  Research  Aircraft  has  a  continuing  requirement  for  more 
accurate,  inertially-based  navigation  data  for  both  track  recovery  and  the  calculation  of  wind  gust 
components.  This  navigational  accuracy  is  necessary,  not  just  during  post-flight  analysis,  but  also  for  real¬ 
time,  in-flight  guidance  and  wind  computation.  Previous  developmental  work  at  the  Flight  Research 
Laboratory  on  advanced  navigation  systems  has  demonstrated  the  benefits  of  a  Kalman  filter  integrated 
navigation  approach  in  order  to  satisfy  the  most  stringent  navigational  requirements.  A  significant  upgrade 
to  the  navigation  sensor  suite  onboard  the  Twin  Otter  in  the  last  two  years  has  resulted  in  the  potential,  via 
Kalman  filtering,  for  generating  very  high  quality  inertial  velocity  and  positional  information  in  real  time, 
together  with  improved  airborne  wind  components. 

The  Kalman  filter  integrated  navigation  design  described  in  this  report  is  based  on  the  optimal 
blending  of  data  from  an  LTN-90-100  strapdown  Inertial  Reference  System  (IRS),  a  Decca  Type  72 
Doppler  velocity  sensing  (DVS)  system  and  an  ARNAV  R-40  airborne  Loran-C  receiver  -  sensors  that  are 
available  on  the  Twin  Otter  at  the  present  time.  In  the  Twin  Otter's  real-time  computing/data  acquisition 
system,  all  three  of  these  navigation  sensors  are  interfaced  to  the  onboard  LSI-1 1/73  microcomputer 
system,  and  a  complete  set  of  navigation  parameters  is  being  recorded.  In  particular,  all  of  the  raw  inertial 
data  parameters  from  the  LTN-90-100  IRS,  required  for  proper  design  of  an  IRS-based  Kalman  filler,  are 
available  with  sufficient  resolution  and  at  a  suitable  digital  sampling  rate. 

A  driving  force  behind  the  decision  to  pursue  this  integrated  navigation  approach  on  the  Twin 
Otter  has  been  the  observation  that  significant  velocity  errors  (and,  eventually,  position  errors)  can  occur 
in  the  LTN-90-100  IRS  over  the  course  of  a  flight,  and  the  observed  error  levels  can  seriously  degrade  the 
accuracy  of  the  wind  calculations.  On  the  other  hand,  the  airborne  Loran-C  positional  data  has  been 
demonstrated  to  be  consistently  more  accurate  than  the  IRS  position  information,  in  the  long  term.  An 
integrated  navigation  system  approach,  using  the  principles  of  Kalman  filtering,  is  shown  to  have  the  ability 
to  use  Loran-C  data  (and,  to  a  lesser  extent,  Doppler  velocity  data)  to  accurately  track  the  dominant  IRS 
errors  (position,  velocity  and  attitude  components),  and  provide  IRS  error  corrections  at  a  rate  appropriate 
for  Twin  Otter  requirements. 
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RESUME 


Le  Twin  Otter  de  I'lRA  charge  des  recherches  atmosph£riques  n£cessite  constamment  une  plus  grande 
precision  des  donnees  de  navigation  par  inertie,  tant  pour  le  r£tablissement  de  la  trajectoire  que  pour  le 
calcul  des  composantes  des  rafales  de  vent.  La  precision  de  la  navigation  est  n£cessaire  non  seulement 
tors  de  I’analyse  consecutive  £  un  vol,  mais  egalement  tors  du  guidage  en  vol  en  temps  reel  et  du  calcul  du 
vent.  Au  Laboratoire  de  recherche  en  vol,  les  systemes  avanc6s  de  navigation  ont  fait  I'objet  de  travaux 
anterieurs  de  mise  au  point  qui  ont  demontre  qu'un  dispositif  de  navigation  integr6e  &  filtre  Kalman 
presentait  des  avantages  permettant  de  satisfaire  aux  exigences  de  navigation  les  plus  severes.  Au  cours 
des  deux  derni6res  ann6es,  une  amelioration  importante  de  I'ensemble  des  detecteurs  installs  £  bord 
du  Twin  Otter  s'est  traduite  par  la  possibility  de  generer,  par  le  biais  d’un  filtre  de  Kalman,  des  donnees  de 
tr£s  haute  qualite  sur  la  vitesse  inertielle  et  £  la  position  en  temps  reel  et  par  ('amelioration  du  calcul  des 
composantes  du  vent  en  cours  de  vol. 

Les  principes  de  navigation  integree  £  filtre  de  Kalman  d£crits  dans  le  present  rapport  sont  fondes 
sur  la  combinaison  optimale  de  donnees  provenant  d'un  syst£me  de  navigation  par  inertie  (SNI)  £ 
composants  li£s  LTN-90-100,  d'un  systfeme  Doppler  de  mesure  de  vitesse  Decca  Type  72  et  d’un 
r£cepteur  Loran-C  de  bord  ARNAV  R-40;  tous  ces  detecteurs  sont  pr£sentement  installs  dans  le  Twin 
Otter.  Le  systdme  de  traitement  en  temps  reel  et  d'acquisition  des  donnees  du  Twin  Otter  est  forme  par 
I'interfagage  de  ces  trois  detecteurs  de  navigation  au  micro-ordinateur  LSI  11/73  de  bord;  le  systdme 
enregistre  done  un  jeu  complet  de  param£tres  de  navigation.  On  dispose  plus  particuli6rement  de  tous 
les  parametres  d'inertie  bruts  du  SNI  LTN-90-100  n£cessaires  £  la  conception  du  type  approprie  de  filtre 
de  Kalman  fonde  sur  SNI,  la  resolution  et  le  taux  d'6chantillonnage  num£rique  etant  suffisants. 

L'une  des  principales  justifications  de  Padoption  de  cette  approche  de  navigation  int6gr£e  sur  le 
Twin  Otter  est  la  possibilite  de  production  pendant  un  vol  d'importantes  erreurs  dans  les  donnees  de 
vitesse  (et,  £ventuellement,  dans  celles  de  position)  du  SNI  LTN-90-100;  les  niveaux  d'erreurs  observes 
peuvent  r£duire  substantiellement  I'exactitude  des  calculs  du  vent.  Par  contre,  on  a  d£montr6  que  les 
donn6es  de  position  du  Loran  C  de  bord  sont  g6n6ralement  plus  exactes  £  tong  terme  que  les  donnees 
de  position  du  SNI.  On  montre  qu'une  approche  de  navigation  int6gr£e  integrant  les  principes  du  filtrage 
de  Kalman  peut  utiliser  les  donnees  Loran-C  (et,  dans  une  moindre  mesure,  les  donnees  de  vitesse 
Doppler)  pour  suivre  avec  precision  les  principales  erreurs  SNI  (composantes  de  position,  vitesse,  et 
attitude)  et  produire  des  corrections  d’erreur  SNI  £  un  debit  satisfaisant  aux  exigences  du  Twin  Otter. 
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A  KALMAN  FILTER  INTEGRATED  NAVIGATION  DESIGN  FOR 
THE  IAR  TWIN  OTTER  ATMOSPHERIC  RESEARCH  AIRCRAFT 


1.0  INTRODUCTION 

The  IAR  Twin  Otter  Atmospheric  Research  Aircraft  has  a  continuing  requirement  for  more 
accurate,  inertially-based  navigation  data  for  both  track  recovery  and  the  calculation  of  wind  gust 
components.  This  navigational  accuracy  is  necessary,  not  just  during  post-flight  analysis,  but  also  for  real¬ 
time,  in-flight  guidance  and  wind  computation.  Previous  developmental  work  at  the  Flight  Research 
Laboratory  (FRL)  on  advanced  navigation  systems  has  demonstrated  the  benefits  of  a  Kalman  filter 
integrated  navigation  approach  in  satisfying  the  most  stringent  navigational  requirements  (Refs.  1 , 2, 3).  A 
significant  upgrade  to  the  navigation  sensor  suite  onboard  the  Twin  Otter  in  the  last  two  years  has  resulted 
in  the  potential,  via  Kalman  filtering,  for  generating  very  high  quality  inertial  velocity  and  positional 
information  in  real  time,  together  with  improved  airborne  wind  components.  Figure  1  shows  the  IAR  Twin 
Otter's  suite  of  navigation,  inertial  sensing  and  air  data  instrumentation. 

The  Kalman  filter  integrated  navigation  design  described  in  this  report  is  based  on  the  optimal 
blending  of  data  from  an  LTN-90-100  strapdown  Inertial  Reference  System  (IRS),  a  Decca  Type  72 
Doppler  velocity  sensing  (DVS)  system  and  an  ARNAV  R-40  airborne  Loran-C  receiver  -  sensors  that  are 
available  on  the  Twin  Otter  at  the  present  time.  In  the  Twin  Otter's  real-time  computing/data  acquisition 
system,  all  three  of  these  navigation  sensors  are  interfaced  to  the  onboard  LSI-1 1/73  microcomputer 
system,  and  a  complete  set  of  navigation  parameters  is  being  recorded.  In  particular,  all  of  the  raw  inertial 
data  parameters  from  the  LTN-90-100  IRS,  required  for  proper  design  of  an  IRS-based  Kalman  filter,  are 
available  with  sufficient  resolution  and  at  a  suitable  digital  sampling  rate. 

A  driving  force  behind  the  decision  to  pursue  this  integrated  navigation  approach  on  the  Twin 
Otter  has  been  the  observation  that  significant  velocity  errors  (and,  eventually,  position  errors)  can  occur 
in  the  LTN-90-100  IRS  over  the  course  of  a  flight,  and  the  observed  error  levels  can  seriously  degrade  the 
accuracy  of  the  wind  calculations.  On  the  other  hand,  the  airborne  Loran-C  positional  data  has  been 
demonstrated  to  be  consistently  more  accurate  than  the  IRS  position  information,  in  the  long  term.  An 
integrated  navigation  system  approach,  using  the  principles  of  Kalman  filtering,  has  the  ability  to  use  Loran- 
C  data  (and,  to  a  lesser  extent,  Doppler  velocity  data)  to  accurately  track  the  dominant  IRS  errors  (in  both 
the  position  and  velocity  components),  and  provide  IRS  error  corrections  at  a  rate  appropriate  for  Twin 
Otter  requirements. 

In  this  report,  Section  2  gives  a  short  description,  in  general  terms,  of  the  mathematical 
nomenclature  and  the  fundamental  equations  used  to  formulate  a  linear,  discrete-time  Kalman  filter.  As 
well,  computationally  efficient  subroutines  for  performing  linear,  discrete-time  Kalman  filtering  are 
described.  Finally,  the  application  of  Kalman  filtering  to  the  general  problem  of  integrated  airborne 
navigation  is  introduced.  Section  3,  along  with  an  associated  set  of  four  appendices,  then  gives  complete 
details  concerning  the  specific  Twin  Otter  IRS/Doppler/Loran-C  Kalman  filter  design  that  is  the  main 
subject  of  this  report. 

A  comprehensive  set  of  software  for  simulating  various  navigation  sensors  has  been  developed  at 
the  FRL  to  accurately  model  realistic  aircraft  flight  trajectory  dynamics  and  error  dynamics  in  sensors  such 
as  medium-accuracy  INS’s,  Doppler  velocity  sensing  (DVS)  systems  and  various  radio  position  fixing 
systems.  Section  4  describes  the  various  software  modules  that  make  up  this  complete  simulation 
package,  as  well  as  the  software  for  performing  the  actual  integrated  navigation  Kalman  filtering.  An 
extensive  series  of  Kalman  filter  runs  using  this  simulation  software  has  demonstrated  the  fundamental 
robustness  and  accuracy  of  the  Kalman  filter  design.  The  simulation  results  indicate  that  there  is  an 
incremental  improvement  to  be  gained  by  using  measurement  prefiltering  on  both  the  Loran-C  and 
Doppler-based  measurements,  but  the  improvement  is  not  considered  significant  enough  to  warrant 
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using  prefiltering  in  a  real-time  application.  Simulation  runs  also  suggest  that,  in  the  presence  of  the 
reasonably  accurate  Loran-C  position-based  measurements,  the  rather  noisy  and  inaccurate  Doppler 
velocity  measurements  contribute  very  little  to  the  estimation  of  the  various  IRS  error  states.  Filter 
robustness  checks  have  confirmed  that  there  is  virtually  no  difference  in  performance  when  using  an  8  Hz 
version  of  the  IRS  digital  data  rather  than  the  usual  1 6  Hz  version,  or  when  updating  the  Kalman  filter  every 
20  seconds  rather  than  every  10  seconds. 

Section  5  describes,  in  some  detail,  Kalman  filter  results  that  are  based  on  two  sets  of  real  flight 
test  data  acquired  specifically  to  assess  filter  performance  under  typical  aircraft  operational  conditions. 
Most  of  the  results  and  conclusions  established  from  the  simulation  studies  have  been  corroborated  from 
an  analysis  of  the  real  flight  data.  In  particular,  it  has  been  confirmed  that  there  is,  indeed,  no  further 
improvement  in  IRS  error  estimation  accuracy  to  be  gained  by  processing  the  Doppler-based 
measurements  along  with  the  Loran-C  measurements.  In  fact,  one  must  be  careful  to  establish  the  proper 
measurement  noise  'weighting'  of  the  Doppler-based  data  so  as  to  avoid  having  it  actually  degrade  filter 
performance  when  used  in  conjunction  with  the  Loran-C  data.  The  best  one  can  achieve  is  an  effective 
ongoing  calibration  of  the  Doppler  velocity  sensing  system  which  would  be  useful  in  a  situation  where 
either  the  IRS  or  Loran-C  sensor  failed  in  flight.  From  the  analysis  of  real  flight  data,  it  has  been  found  that 
under  certain  circumstances,  the  use  of  Visual  On-Tops  (VOT's)  as  extra  measurements  can  cause  an 
undesirable  transient  in  several  of  the  Kalman  filter  error  states.  The  procedures  developed  for  using 
VOT's  as  measurements  must  be  checked  out  very  carefully,  for  any  given  Kalman  filter  configuration,  to 
ensure  that  basic  filtering  integrity  is  maintained. 

General  conclusions  and  future  work  are  discussed  in  Section  6.  The  application  of  Kalman 
filtering  techniques  to  the  analysis  of  Twin  Otter  flight  data  is  an  evolving  process  with  several  hardware 
and  software  changes  planned  for  the  near  future  as  outlined  in  Section  6. 


2.0  THE  KALMAN  FILTER  EQUATIONS  FOR  INTEGRATED  NAVIGATION 

2.1  The  Linear  Discrete-Time  Kalman  Filter 

For  completeness,  and  to  establish  the  nomenclature  that  will  be  used  throughout  the  report,  the 
general  form  of  the  linear,  discrete-time  version  of  the  Kalman  filter  equations  will  be  outlined  here. 
Readers  are  directed  to  Refs.  4,  5  and  6  for  a  thorough  discussion  of  the  fundamental  theory  of  Kalman 
filtering. 


Assume  that  a  physical  system  has  an  equivalent  n*h-order,  discrete-time  dynamic  model  of  the 

form 


x(k+1)  =  <D(k,k+1)  x(k)  +  G(k)  u(k) 


(1) 


where  x(k)  is  the  n'h-order  system  state  vector  evaluated  at  discrete  time  t|<;  4>(k,k+l)  is  the  n  x  n  state 
transition  matrix  over  the  time  interval  tk  ~>  t^+i ;  G(k)  is  the  n  x  r  plant  noise  gain  matrix  at  tk ;  and  u(k)  is  the 
reorder  vector  of  zero-mean,  white,  Gaussian  (ZMWG)  discrete  plant  noise  processes  having  covariance 
matrix  Q(k)  at  tk- 

For  the  above  dynamic  system,  let  a  discrete-time,  m^-order  measurement  process  exist  in  the 

form 


z(k+l) 


H(k+1)  x(k+1)  +  v(k+1) 


(2) 
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where  z(k+1)  is  the  measurement  vector  at  time  t^+i,  H(k+1)  is  the  m  x  n  observation  matrix  and  v(k+l)  is 
the  mfri-order  measurement  noise  vector  having  covariance  matrix  R(k+1)  at  t|<+i .  Assume  that  noise 
vectors  u  and  v  are  statistically  independent  (i.e.  the  components  of  u  are  uncorrelated  with  the 
components  of  v);  and  assume  also  that  x(0)  is  independent  of  both  u  and  v.  Define  the  following  vector 
and  matrix  variables: 

x'(k+l ),  P'(k+1 )  -  time  update  of  the  state  vector  and  its  covariance  at  tk+1 

(i.e.  just  prior  to  a  measurement  update  at  time  tk+i)- 

x(k+l ),  P(k+1 )  -  optimal  state  estimate  and  its  associated  covariance  at  t^+i 

(i.e.  just  after  a  measurement  update  at  time  t^+i). 

xo,  Po  -  initial  conditions  on  the  state  vector  and  its  covariance  (i.e.  x(0)  =  Xo,  P(0)  =  Po). 

Under  the  foregoing  definitions  and  assumptions,  it  can  be  shown  that  the  optimal  estimate  of  the 
state  vector  at  time  tk+i  (i.e.  x(k+1))  and  its  associated  error  covariance  (i.e.  P(k+1 ))  can  be  computed  from 
the  following  set  of  five  recursion  equations  that  form  the  heart  of  discrete-time  Kalman  filtering: 


x'(k+1)  =  <&(k,k+1 )  x(k)  state  time  update 

P'(k+1)  =  fl>(k,k+1)  P(k)  <DT(k,k+1)  +  G(k)  Q(k)  GT(k)  error  covariance  time  update 

K(k+1)  =  P'(k+1)H(k+1)T  [H(k+1)P'(k+1)HT(k+1)  +  R(k+1)]-i  Kalmangain 

x(k+1)  =  x'(k+1)  +  K(k+1)[z(k+1)  -  H(k+1)  x'(k+1)]  state  measurement  update 

P(k+1)  =  P'(k+1)  -  K(k+1)  H(k+1)  P'(k+1)  error  covariance  measurement  update 

(3) 

with  initialization  of  this  recursive  procedure  provided  by  a  priori  knowledge  of  x0  and  Po  It  should  be 
noted  that,  when  an  error  state  Kalman  filter  is  being  used  (see  subsection  2.3),  the  state  vector  in  Eqns. 
3  is  denoted  as  8x. 


2.2  Software  for  Performing  Discrete-Time  Kalman  Filtering 

Numerically  efficient  software  exists  for  performing  the  foregoing  set  of  five  recursion  equations, 
especially  if  certain  simplifying  assumptions  can  be  made.  This  author  uses  a  set  of  four  relatively  short 
subroutines  to  perform  the  above  veotor/matrix  calculations  in  a  very  efficient  fashion.  They  are  based  on 
Bierman’s  UDUT  factorization  algorithms,  which  are  generally  acknowledged  to  be  the  most  numerically 
stable,  computationally  efficient  ones  to  use,  especially  for  real-time  applications.  In  order  to  avoid  the 
matrix  inversion  that  would  normally  be  required  in  the  Kalman  gain  equation  of  Eqns.  3,  the  assumption  is 
made  that  the  elements  of  z,  the  measurement  vector,  are  statistically  independent  of  each  other.  This 
assumption  is  certainly  valid  for  the  kinds  of  measurements  that  are  described  in  subsection  3.3.  Under 
this  assumption,  the  covariance  matrix  R  is  diagonal  and  it  can  be  shown  (Ref.  7)  that  the  measurement 
vector  z  can  then  be  processed  in  a  one-at-a-time  fashion.  In  essence,  the  Kalman  gain  vector/matrix 
equation  is  converted  to  a  sequence  of  much  simpler  scalar/vector  operations.  The  four  subroutines  that 
make  up  the  UDUT  factorization  version  of  Kalman  filtering  are  described  as  follows: 
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I)  SUBROUTINE  FACTOR(N,NN,P,U)  -  a  subroutine  to  factor  the  covariance  matrix  P  into 
matrices  U  and  D,  where  P  =  UDUT.  U  is  upper  triangular  with  unit  diagonal,  and  D  is  a  diagonal  matrix. 
The  inputs  to  the  subroutine  are  as  follows: 

N  -  dimension  of  the  state  vector  x; 

NN  -  state  dimension  augmented  by  plant  noise  dimension 
i.e.  NN  =  N  +  NW,  where  NW  is  size  of  u; 

P(N,N)  -  covariance  matrix  to  be  factored. 


The  outputs  from  the  subroutine  are  the  following: 

U(N,NN)  -  the  matrix  of  the  factors  of  P,  with  D  on  the  diagonal  and  extra  storage  for  matrix  G; 
P(N,N)  -  modified  covariance  matrix,  with  upper  triangular  portion  overwritten. 

II)  SUBROUTINE  UNFACTfN.NN.P^)  -  this  subroutine  reconstructs  the  P  matrix  from  its 
UDUT  decomposition.  The  inputs  to  the  subroutine  are  the  following: 

N  -  dimension  of  the  state  vector  x; 

NN  •  state  dimension  augmented  by  plant  noise  dimension; 

U(N,NN)  -  the  matrix  of  the  factors  of  P,  with  0  on  the  diagonal  and  extra  storage  for  G. 

The  outputs  from  the  subroutine  are: 

P(N,N)  -  reconstructed  covariance  matrix. 

iii)  SUBROUTINE  MWGSUD  (N.NW.NN.U.G.PHI.Q.X.A.V.D)  -  this  subroutine  performs 
Bierman's  modified  weighted  Gram-Schmidt  (MWGS)  time  update  of  state  estimate  x  and  its  factored 
covariance  matrix  U  (i.e.  the  first  two  equations  in  the  five-equation  set  of  Eqns.  3).  The  inputs  to  this 
subroutine  are: 

N  -  dimension  of  the  state  vector  x; 

NW  -  dimension  of  the  plant  noise  vector  u; 

NN  -  augmented  state  dimension,  where  NN  =  N  +  NW; 

U(N,NN)  -  a  matrix  that  contains  the  factored  covariance  U(N,N)  in  the  first  N  columns  and 
G(N,NW)  in  the  remaining  NW  columns; 

G(N,NW)  -  plant  (process)  noise  gain  matrix  (G  in  Eqn.  1); 

PHI(N,N)  •  linear  state  transition  matrix  ($  in  Eqn.  1); 

Q(NW)  -  vector  of  plant  noise  covariance  diagonal  elements; 

X(N)  -  optimal  state  estimate  from  previous  Kalman  filter  iteration. 
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The  outputs  from  this  subroutine  are  the  following: 

X(N)  -  time  updated  state  estimate  (first  of  Eqns.  3); 

U(N,NN)  -  time  updated  covariance  factors  (second  of  Eqns.  3); 

A(NN),  V(NN),  D(NN)  -  double  precision  working  vectors. 

Iv)  SUBROUTINE  UDUPDCfN.NN.X.U.ZZ.RR.H.K.RSD.ALPHA, ITST, IPASS, A,B)  - 
this  subroutine  performs  Bierman's  UDUT  measurement  update  for  a  scalar  measurement  z  (i.e.  the  last 
three  of  Eqns.  3).  Prior  to  updating  U  and  x,  a  tolerance  test  can  be  performed  on  the  measurement 
residual.  The  measurement  residuals  sequence  (also  called  the  innovations)  and  an  estimate  of  its 
variance  are  computed  as  well.  If  the  residual  lies  outside  of  the  1  -sigma  boundaries  defined  by  the 
tolerance,  the  subroutine  returns  without  updating  U  and  x.  The  inputs  to  this  subroutine  are: 

N  -  dimension  of  state  vector  x; 

NN  -  state  dimension  augmented  by  plant  noise  dimension; 

X(N)  -  state  vector  from  a  time  update  or  intermediate  measurement  update; 

U(N,NN)  -  factored  covariance  matrix  from  a  time  update  or  intermediate  measurement  update; 

ZZ  -  scalar  measurement  z  to  be  processed; 

RR  -  noise  variance  associated  with  the  measurement; 

H(N)  -  geometry  vector  for  the  measurement  z  (i.e.  appropriate  row  of  observation  matrix  H); 

ITST  -  integer  multiplier  to  set  tolerance  for  residuals  monitoring.  If  ITST  is  set  negative,  then 
residuals  monitoring  is  omitted. 

The  outputs  from  the  subroutine  are  as  follows; 

X(N)  -  updated  state  vector  based  on  the  measurement; 

U(N,NN)  -  updated  covariance  matrix  in  factorized  form; 

K(N)  -  Kalman  gain  vector  associated  with  the  measurement  (i.e.  appropriate  column  of  Kalman 
gain  matrix  K  in  the  third  one  of  Eqns.  3); 

RSD  -  latest  value  of  measurement  residuals  (innovation); 

ALPHA  -  latest  estimate  of  the  variance  of  the  residuals  sequence; 

IPASS  -  integer  flag  to  indicate  if  measurement  is  used  (IPASS  =  1)  or  not  used  (IPASS  =  0) 
when  residuals  monitoring  is  invoked; 


A(N),  B(N)  -  working  vectors. 
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2.3  Integrated  Navigation  Using  Kalman  Filtering 

There  are  three  generally  accepted  configurations  for  the  use  of  Kalman  filtering  to  blend,  in  an 
optimal  fashion,  navigation  data  from  various  airborne  transducers.  These  three  configurations  (Refs.  3, 8) 
are  known  as  the  i)  total  state  representation,  the  ii)  error  state  feedforward  representation  and  the  iii)  error 
state  feedback  representation.  In  all  three  cases,  navigation  data  from  a  dead  reckoning  (DR)  navigator 
such  as  inertial  or  Doppler/heading  reference  are  blended  together  with  redundant  navigation  data  from 
one  or  more  external  navaids  such  as  Omega/VLF,  airborne  Loran-C,  Navstar  GPS  or  VOR/DME.  The 
optimal  blending  is  accomplished  in  real  time  through  the  use  of  a  Kalman  filter  algorithm  that  differs 
somewhat  depending  upon  the  configuration  being  used. 

In  this  report,  the  error  state  feedforward  configuration  will  be  used  exclusively.  Figure  2  shows  a 
simplified  block  diagram  representation  of  this  configuration.  In  the  error  state  feedforward 
implementation,  the  inertial  navigation  system  (INS)  or  Doppler/heading  dead  reckoning  (DR)  system  is 
treated  as  a  'black  box'  that  outputs  position,  velocity  and  (possibly)  attitude  data  corrupted  with  errors. 
Measurements,  in  the  form  of  position,  velocity  or  range  differences  between  DR  predictions  and  what  is 
actually  being  sensed  by  an  external  navaid,  are  sent  to  the  Kalman  filter  along  with  fundamental  data 
generated  by  the  DR  device.  The  Kalman  filter  then  estimates  the  dominant,  low  frequency  error  states  of 
the  DR  system  and  the  external  navaids  in  real  time  (i.e.  as  they  are  occurring).  The  DR  errors  are  nulled 
out,  in  a  feedforward  sense,  to  produce  best  estimates  of  navigation  variables  such  as  3-D  position, 
velocity  and  the  Euler  angles  (i.e.  0,  <t>  and  y). 

It  is  worthwhile  to  list  the  salient  features  of  the  error  state  feedforward  approach  to  Kalman  filter 
integrated  navigation: 

Any  potentially  nonlinear  navigation  calculations  in  the  DR  'black  box'  are  performed  separately 
from  the  Kalman  filter  recursive  update  calculations. 

Any  required  nonlinear  pre-processing  of  external  navaid  data  also  occurs  separately  from  the 
Kalman  filter  equations.  Position,  velocity  and/or  range  differences  between  DR  predictions  and  navaid 
outputs  are  calculated  and  then  sent  to  the  Kalman  filter  as  measurement  data. 

Under  the  foregoing  conditions,  the  Kalman  filter  can  be  effectively  formulated  as  a  sub-optimal, 
reduced-order,  linear  Kalman  filter  (Ref.  8)  which  will  estimate  the  dominant  navigation  error  states,  as 
opposed  to  the  actual  physical  states  of  the  navigation  system.  Also,  by  using  the  error  state  formulation, 
the  use  of  a  linear  Kalman  filter  design  is  more  likely  to  be  a  valid  assumption. 

The  error  states  and  the  measurements  used  in  the  Kalman  filter  are  all  relatively  small  in 
magnitude,  so  numerical  round-off  error  is  usually  not  a  problem  even  when  cycling  through  the  recursion 
equations  (i.e.  Eqns.  3)  many  times. 

The  dominant  error  states  of  any  navigation  sensor  vary  slowly  in  time  relative  to  the  actual  sample 
rate  of  the  navigation  data.  Hence,  discrete-time  Kalman  filtering  can  take  place  at  a  slow  update  rate;  but 
the  error  states  so  estimated  can  be  applied  to  the  navigation  data  at  the  much  faster  navigation  data  rate. 

With  a  discrete-time,  linear  Kalman  filter  design,  efficient  'off-the-shelf'  FORTRAN  software,  such 
as  Bierman's  UDUT  factorization  algorithms  described  in  subsection  2.2,  can  be  used  (Ref.  7). 
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3.0  IRS/DOPPLER/LORAN-C  KALMAN  FILTER  NAVIGATOR 

3.1  Navigation  Sensor  Descriptions 

At  the  present  time,  the  Twin  Otter's  suite  of  navigation  sensors  (Fig.  1)  available  for  optimal 
blending  is  as  follows: 

I)  Litton  LTN-90-100  inertial  Reference  System  -  this  modern,  three-axis  strapdown  inertial 
system  (called  an  IRS  -  Inertial  Reference  System  -  by  Litton)  utilizes  ultra-reliable  ring  laser  gyro 
technology  and  a  digital  data  bus  (ARINC  429  standard)  to  provide  a  complete  set  of  inertial  parameters  in 
digital  form,  many  at  a  64  Hz  datarate.  Table  3.1  gives  a  list  of  the  16  inertial  parameters  recorded  from  the 
digital  data  bus  (and  required  for  the  IRS-based  Kalman  filter  design)  together  with  the  update  rate,  units, 
dynamic  range,  resolution  and  positive  sense  for  each  parameter.  The  specific  interface  design  for  the 
Twin  Otter  onboard  data  acquisition  system  is  such  that  all  16  IRS  parameters  are  sampled  at  16  Hz.  Table 

3.2  lists  the  fundamental  specifications  of  the  system's  accelerometers  and  gyroscopes  -  these  statistical 
specifications  will  be  of  importance  in  the  Kalman  filter  design.  This  particular  IRS  requires  barometric 
altimeter  data,  in  ARINC  429  format,  as  an  input  in  order  to  stabilize  the  vertical  channel  via  an  internal,  third- 
order,  fixed-gain  digital  filter.  Typical  accuracy  for  this  type  of  IRS  would  be  1  nm/hr  in  the  horizontal 
position  components,  5  knots  in  the  horizontal  velocity  components,  0.05  deg  for  the  pitch  and  roll 
attitude  components  and  0.2  deg  for  the  heading.  The  reader  is  referred  to  Litton's  Technical  Description 
Manual  (Ref.  9)  for  the  LTN-90-100  IRS  for  further  technical  details  about  this  IRS. 

II)  Decca  Doppler  Radar  Type  72  -  this  3-beam  Janus  Doppler  radar  system  measures  the 
three  body-axis  (i.e.  strapdown)  components  of  aircraft  velocity,  namely  U,  V  and  W.  The  specially- 
designed  digital  interface  for  this  unit  is  such  that  each  Doppler  velocity  component  is  averaged  over  a  1/2 
second  interval  and  then  sampled,  thus  yielding  an  effective  sample  rate  of  2  Hz.  Table  3.3  shows  some 
of  the  important  characteristics  of  this  Doppler  radar  system  as  it  is  interfaced  to  the  Twin  Otter  onboard 
data  acquisition  system,  including  the  range,  accuracy  and  resolution  of  each  of  the  velocity  components. 
For  further  information  on  the  technical  details  of  this  Doppler  system,  the  reader  should  consult  Ref.  10. 

III)  ARNAV  R-40-AVA-1000A  Loran-C  Receiver  -  this  airborne  Loran-C  receiver  provides 
digitized  geographical  latitude  and  longitude  data,  at  a  nominal  1  Hz  update  rate,  via  a  standard  RS-232C 
serial  output  port.  Within  the  normal  coverage  area  of  one  of  the  Loran-C  station  chains,  expected 
accuracy  is  0.2  nm  or  better  in  each  of  latitude  and  longitude.  The  resolution  of  the  Loran-C  latitude  and 
longitude  data,  as  digitized  onboard  the  Twin  Otter,  is  0.01  arc  minutes.  The  Operation  Manual  for  the  R- 
40  Loran-C  receiver  (Ref.  11)  provides  technical  information  on  the  Loran-C  system  as  a  whole,  as  well  as 
on  the  R-40  receiver  itself.  As  an  indication  of  expected  airborne  Loran-C  performance,  Ref.  12  shows 
results  of  flight  tests  that  took  place  during  the  certification  of  an  airborne  Loran-C  navigation  system. 


3.2  Overview  of  the  Kalman  Filter  Navigator 

A  simple  block  diagram  representation  of  the  proposed  IRS/Doppler/Loran-C  Kalman  filter 
navigator  is  shown  in  Fig.  3.  This  figure  shows  a  specific  implementation  of  the  error  state  feedforward 
version  of  a  Kalman  filter  integrated  navigation  scheme.  In  this  particular  design,  the  LTN-90-100  IRS  is 
the  dead  reckoning  (DR)  system  that  outputs  error-corrupted  position,  velocity  and  attitude  information  at 
a  16  Hz  datarate.  Also  shown  is  the  air  data  source  of  altitude  information  that  is  required  by  the  IRS  for  its 
third-order  baro  damping  loop  (see  Appendix  A  for  loop  details).  Two  other  navaids,  the  airborne  Loran-C 
receiver  and  the  Doppler  radar  system,  supply  redundant  navigation  information  that  is  used  to  form 
measurements  to  be  processed  by  the  error  state  Kalman  filter.  In  the  filter  design  shown  here,  the 
differences  in  position  (both  latitude  and  longitude  components)  between  the  IRS  and  the  Loran-C 
receiver  are  used  as  measurements  for  the  Kalman  filter.  This  particular  measurement  type  could  be 
processed  by  the  Kalman  filter  as  frequently  as  1  Hz,  since  that  is  the  rate  at  which  Loran-C  data  is 
sampled. 
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TABLE  3.1 

DIGITAL  DATA  PARAMETERS  AVAILABLE  FROM  LTN-90-100  IRS 


PARAMETER 

RATE 

(HZ) 

UNITS 

RANGE(+/-) 

RESOLUTION 

VVE  SENSE 

1. 

Latitude 

8 

Degs 

180 

.000172 

North  Fr.  0  0 

2. 

Longitude 

8 

Degs 

180 

.000172 

East  Fr.  0  0 

3. 

True  Heading 

32 

Degs 

180 

.005493 

CW  Fr.  North 

4. 

Pitch  Angle 

64 

Degs 

180 

.005493 

Up 

5. 

Roll  Angle 

64 

Degs 

180 

.005493 

Rgt  Wing  Down 

6. 

Body  Ptch  Rate 

64 

Deg/S 

128 

.003906 

Up 

7. 

Body  Roll  Rate 

64 

Deg/S 

128 

.003906 

Rgt  Wing  Down 

8. 

Body  Yaw  Rate 

64 

Deg/S 

128 

.003906 

Nose  Right 

9. 

Body  Long  Acc 

64 

G 

4 

.000122 

Forward 

10. 

Body  Lat  Accel 

64 

G 

4 

.000122 

Right 

11. 

Body  Norm  Acc 

64 

G 

4 

.000122 

Up 

12. 

Inertial  Alt'd 

32 

Feet 

131,072 

0.125 

Up 

13. 

Vertical  Accel 

64 

G 

4 

.000122 

Up 

14. 

Inert  Vert  Spd 

32 

Ft/M  in 

32,768 

1.00 

Up 

15. 

N-S  Velocity 

16 

Knots 

4096 

0.125 

North 

16. 

E-W  Velocity 

16 

Knots 

4096 

0.125 

East 

TABLE  3.2 

SPECIFICATIONS  FOR  IRS  ACCELEROMETERS  AND  GYROSCOPES 


A-4  ACCELEROMETERS; 

Scale  Factor  Repeatability 
Scale  Factor  Nonlinearity 
Bias  Repeatability 
Bias  Short  Term  Stability 
Alignment  Stability 
Maximum  Acceleration 


0.005%  (1 -sigma, 1  year) 
10  pg/g  2(1  -sigma) 

50  pg(1  -sigma,  1  year) 

5  ng  (1 -sigma) 

5  arc  sec  (1 -sigma) 

25  g 


LG-8Q28B  GYROS; 

Scale  Factor  Repeatability 
Bias  Repeatability 
Random  Drift 
Maximum  Rate 


0.0005%  (1 -sigma) 

0.01  deg/hr  (1 -sigma) 
0.003  deg/Vhr  (1 -sigma) 
400  deg/s 
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TABLE  3.3 

DIGITAL  DATA  PARAMETERS  AVAILABLE  FROM  DECCA  DOPPLER  SYSTEM 


COMPONENT 

UNITS 

RANGE 

RESOLUTION 

ACCURACY 

+'V  SENSE 

Long  Vel  (U) 

Knots 

497.3 

0.0194 

2.0 

Forward 

Lat  Vel  (V) 

Knots 

248.6 

0.0097 

2.0 

Right 

Norm  Vel  (W) 

Knots 

124.3 

0.0039 

2.0 

Down 

The  second  type  of  measurement  data  available  for  processing  is  derived  from  the  Doppler 
strapdown  velocity  components.  As  can  be  seen  in  Fig.  3,  Doppler  U.V.W  velocity  components  are 
differenced  with  their  IRS  counterparts,  and  these  strapdown  (i.e.  body  axis)  velocity  differences  are  then 
processed  by  the  Kalman  filter.  The  data  rate  for  velocity  measurements  is  the  fundamental  Doppler  rate 
of  2  Hz.  In  the  error  state  Kalman  filter  design,  dominant  sources  of  error  in  the  IRS,  Loran-C  receiver  and 
Doppler  radar  system  are  modelled  and  estimated  in  an  on-line  fashion  as  the  measurement  data  are 
processed.  A  simple  corrector  algorithm  is  then  used  to  transform  tilt  angle  error  estimates  into  Euler 
angle  errors,  and  the  16  Hz  data  stream  of  raw  IRS  data  (position,  velocity  and  attitude  components)  is 
corrected,  in  real  time,  using  the  slowly  varying  error  estimates  generated  by  the  Kalman  filter.  The  basic 
update  interval  for  the  Kalman  filter  processing  (i.e.  executing  Eqns.  3  for  a  new  set  of  measurement  data) 
is  set  at  ten  seconds  -  more  than  adequate  for  'tracking'  the  expected  sources  of  error  in  the  various 
navigation  systems.  Not  only  is  the  IRS  corrected  in  a  feedforward  sense  -  the  Loran-C  receiver  and 
Doppler  radar  system  are  corrected  for  their  dominant  low  frequency  errors  as  well.  This  on-line  calibration 
of  the  redundant  navaids,  as  well  as  the  IRS,  allows  for  the  possibility  of  a  'graceful  degradation*  in 
navigational  accuracy  should  the  main  IRS-based  DR  system  fail  in  flight. 


3.3  Details  of  the  Kalman  Filter  Design 
3.3.1  Error  States  Chosen 

For  the  complete  IRS/Doppler/Loran-C  Kalman  filter,  there  are  a  total  of  24  error  states  modelled 
(Appendix  B  shows  details  of  the  error  modelling  for  the  IRS,  Doppler  and  Loran-C  systems).  The  error 
states  are  divided  into  two  groups  -  i)  10  system  error  states  that  relate  to  the  basic  baro-damped  IRS  (i.e. 
errors  in  inertial  position,  velocity,  attitude  components  and  vertical  loop  acceleration  error);  and  ii)  14  first- 
order  Markov  error  states  that  correspond  to  the  slowly  varying,  bias-like  errors  assumed  to  exist  in  the 
inertial  sensors  (i.e.  bias  errors  in  accelerometers,  gyros,  altimeter)  and  redundant  navaids  (i.e.  Loran-C 
offsets,  Doppler  scale  factor  and  boresight  errors,  sea  currents).  Table  3.4  lists  all  of  the  system  error 
states,  together  with  typical  values  of  the  1  -sigma  levels  that  would  be  used  for  the  initial  conditions,  Po, 
when  running  the  Kalman  filter.  Table  3.5  is  a  list  of  the  14  Markov  error  states,  again  showing  typical  1- 
sigma,  initial  condition  values  as  well  as  the  nominal  values  of  correlation  times  that  would  be  assumed  for 
the  associated  first-order  Markov  error  processes.  It  should  be  noted  that  the  last  two  error  states  in  Table 
3.5,  the  two  sea  bias  components,  would  only  be  included  in  the  Kalman  filter  in  the  case  of  an  overwater 
flight.  Note,  also,  that  in  these  tables  the  engineering  units  shown  correspond  to  those  being  used 
internally  by  the  Kalman  filter  and  are  Qfll  metric  units,  as  it  turns  out. 
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3.3.2  Plant  Dynamics  (<&,  G,  Q) 

The  full  24  x  24  discrete  slate  transition  matrix  for  the  Kalman  filter  design  is  derived  from  the 
continuous-time  version  of  the  error  state  equations,  the  details  of  which  have  been  developed  in 
Appendix  B.  The  general  relationship  between  the  continuous-time  matrices,  F(t)  and  G(t),  (shown  in 
Eqn.  B1  of  Appendix  B)  and  their  discrete-equivalent  counterparts,  <t>(k,k+1)  and  G(k),  can  be  expressed 
as  follows: 

tk+1 

Oj  j(k,k+1)  =  /  Fj  j(t)dt ;  i,j  =  1 . 24;  i  *  j 

tk 

tk+1 

0|  i(k,k+1)  =  1  +  1  Fj  j(t)dt ;  i  =  1.....10 

tk 

Ojj(k,k+1)  =  exp(-AT/tj);  i  =  11 . 24 


tk+1 

G  i,j(k)  =  (1/AT)  I  G  ij(t)dt ;  i,j  =  1 . 24  (4) 

tk 


where  AT  =  tk+1  -  tk  is  the  Kalman  filter  update  interval  of  ten  seconds,  the  tj’s  are  the  Markov  error  state 
correlation  times  and  the  integrations  are  performed  numerically  using  a  simple  trapezoidal  integration 
algorithm. 


TABLE  3.4 

IRS  SYSTEM  ERROR  STATES  AND  STATISTICS 


ERROR 

STATE 

DESCRIPTION 

INITIAL 

RMS  VALUE 

UNITS 

1.  8L 

Latitude  Error 

2.91  x  10-5 

Rad 

2.  8X 

Longitude  Error 

4.11  x  10-5 

Rad 

3.  8h 

Altitude  Error 

30 

Feet 

4.  8vn 

North  Velocity  Error 

0.33 

Ft/Sec 

5.  8ve 

East  Velocity  Error 

0.33 

Ft/Sec 

6.  8vz 

Vertical  Velocity  Error 

0.33 

Ft/Sec 

7.  EN 

N  Axis  Attitude  Error 

5.0  x  10-5 

Rad 

8.  eg 

E  Axis  Attitude  Error 

5.0  x  10-5 

Rad 

9 .  £z 

z  Axis  Attitude  Error 

1.5  x  10-3 

Rad 

10. 8a 

V  Accel  Correction 

0.25  x  10-1 

Ft/S2 
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TABLE  3.5 


MARKOV  ERROR 

STATES  AND 

STATISTICS 

ERROR 

DESCRIPTION 

INITIAL 

UNITS 

CORREL 

STATE 

RMS  VALUE 

TIME 

1 1 .  Bcox 

X  Axis  Gyro  Bias 

9.0  x  10-8 

Rad/S 

7,500  S 

12.  B(Oy 

Y  Axis  Gyro  Bias 

9.0  x  10-8 

Rad/S 

7,500  S 

13.  B(oz 

Z  Axis  Gyro  Bias 

9.0  x  10-8 

Rad/S 

7,500  S 

14.  Bax 

X  Axis  Accel  Bias 

2.0  x  10-3 

Ft/S2 

15,000  S 

15.  Bay 

Y  Axis  Accel  Bias 

2.0  x  10-3 

Ft/S2 

15,000  S 

16.  Baz 

Z  Axis  Accel  Bias 

2.0  x  10-3 

Ft/S2 

15,000  S 

17.  Bhb 

Baro  Altimeter  Bias 

600 

Feet 

5,000  S 

18.  Blat 

Loran  Latitude  Bias 

8.73  x  10-5 

Rad 

15,000  S 

19.  Blng 

Loran  Longitude  Bias 

1.23x10-4 

Rad 

15,000  S 

20.  SF(j 

Doppler  SF  Error 

2% 

— 

15,000  S 

21.  By 

V  Dop  Boresight  Error 

0.02 

Rad 

15,000  S 

22.  Bw 

W  Dop  Boresight  Error 

0.02 

Rad 

15,000  S 

23.  SBn 

N  Sea  Bias  Component 

4.5 

Ft/S 

900  S 

24.  SBe 

E  Sea  Bias  Component 

4.5 

Ft/S 

900  S 

A  complete  list  of  the  discrete-equivalent  ZMWG  noise  components  included  in  the  Kalman  filter 
design  is  shown  in  Table  3.6,  along  with  typical  values  for  the  associated  standard  deviations.  Note  that 
the  first  ten  components  in  Table  3.6  correspond  to  integrated  random  noise  effects  over  the  Kalman  filter 
update  interval  (i.e.  AT)  of  ten  seconds.  The  rest  of  the  noise  components  in  Table  3.6  are  discrete- 
equivalent  versions  of  the  so-called  Markov  plant  noises  (Ref.  8)  associated  with  error  states  #1 1  to  #24  in 
Table  3.5.  Their  standard  deviations  are  determined  from  the  filter  update  interval  (AT),  plus  the 
information  given  in  Table  3.5  for  the  correlation  time  and  initial  RMS  level  of  the  associated  Markov  error 
state.  For  simplicity  in  the  filter  design,  all  plant  noise  components  are  assumed  to  have  constant 
variances  that  are  not  affected  by  aircraft  manoeuvring.  This  results  in  a  plant  noise  covariance  matrix  Q 
that  is  diagonal  and  has  all  components  constant  with  time. 

3.3.3  Measurement  Process  (H,  R) 

Recall  that  there  are  two  types  of  measurements  to  be  processed  by  the  Kalman  filter, 
namely  i)  measurements  based  on  Loran-C  position  data,  and  ii)  measurements  based  on  Doppler 
velocity  data.  Details  of  the  mathematical  developments  for  each  type  of  measurement  are  shown  in 
Appendix  C,  while  only  the  important  end  results  are  given  here. 

Measurements  Based  On  Loran-C  - 

These  measurements,  taken  at  discrete  times  tx,  will  consist  of  the  two  simultaneous, 
independent  difference  quantities,  zi(k)  and  z2(k),  where 


Zi(k)  =  LAT|Rs(tk)  -  LATuofl(t|0 

z2(k)  =  LMG|Rs(t|0  -  LNGLOr(!k)  (5) 
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TABLE  3.6 

PLANT  NOISE  COMPONENTS  AND  STATISTICS 


COMPONENT 

DESCRIPTION 

RMS  VALUE 

UNITS 

1. 

U(0X 

X  Gyro  Random  Drift 

4.0  x  10-8 

Rad 

2. 

UG)y 

Y  Gyro  Random  Drift 

4.0  x  10-8 

Rad 

3. 

utoz 

Z  Gyro  Random  Drift 

4.0  x  10-8 

Rad 

4. 

uax 

X  Accel  Random  Noise 

1.0  x  10-3 

Ft/S 

5. 

uay 

Y  Accel  Random  Noise 

1.0  x  10-3 

Ft/S 

6. 

uaz 

Z  Accel  Random  Noise 

1.0  x  10-3 

Ft/S 

7. 

ugN 

N  Random  Gravity 

1.0  x  10-3 

Ft/S 

8. 

ugg 

E  Random  Gravity 

1.0  X  10-3 

Ft/S 

9. 

ugz 

Z  Random  Gravity 

1.0  X  10-3 

Ft/S 

10. 

uhb 

Altimeter  Random  Noise 

31.62 

Ft-S 

11. 

uBo)x 

X  Gyro  Markov  Plant  Noise 

0.46  x  10-8 

Rad/S 

12. 

UBtOy 

Y  Gyro  Markov  Plant  Noise 

0.46  x  10  8 

Rad/S 

13. 

uBwz 

Z  Gyro  Markov  Plant  Noise 

0.46  x  10-8 

Rad/S 

14. 

uBax 

X  Acc  Markov  Plant  Noise 

7.3  x  10  5 

Ft/S2 

15. 

uBay 

Y  Acc  Markov  Plant  Noise 

7.3  x  10-5 

Ft/S2 

16. 

uBaz 

Z  Acc  Markov  Plant  Noise 

1.46  X  10  4 

Ft/S2 

17. 

uBhb 

Altim  Markov  Plant  Noise 

37.91 

Feet 

18. 

uBlt 

Loran  Lat  Markov  Plant  Noise 

0.32  x  10  5 

Rad 

19. 

uBLg 

Loran  Lng  Markov  Plant  Noise 

0.45  x  10-5 

Rad 

20. 

uSFu 

Dop  SF  Markov  Plant  Noise 

0.73  x  10-3 

— 

21. 

uBv 

V  Boresight  Markov  Noise 

0.73  x  10-3 

Rad 

22. 

uByv 

W  Boresight  Markov  Noise 

0.73  x  10-3 

Rad 

23. 

uSBn 

N  Sea  Bias  Markov  Noise 

0.667 

Ft/S 

24. 

uSBg 

E  Sea  Bias  Markov  Noise 

0.667 

Ft/S 

with  LATirs.  LNGirs  the  IRS  estimates  of  geographical  position  and  LATlor.  LNGlor  the  Loran-C 
estimates  of  position.  The  24-element  rows  of  the  H  observation  matrix  that  are  associated  with  these  two 
measurements  (i.e.  the  first  two  rows  of  H)  are  given  very  simply  as  (see  Appendix  C  for  details): 


Hi  =  (1  0  0  0.. 

. . .  -1  0  0  0  0  0  0)  ; 

’-V  is  18*1  element 

H2  =  (010  0... 

..0-1  00000); 

'-1'  is  19*1  element 

(6) 

The  noise  variances  assumed  for  these  two  measurements,  r,  and  r2,  are  the  first  two  diagonal 
elements  of  the  overall  measurement  noise  covariance  matrix  R.  Nominal  values  of  n  =  (0.07  arc  min)2  and 
xz  =  (0.1  arc  min)  2  are  assigned  to  these  parameters,  based  on  the  observed  performance  of  the  Loran-C 
receiver  onboard  the  IAR  Twin  Otter.  The  Loran-based  measurement  data  are  available  at  a  1  Hz  rate;  but 
the  normal  Kalman  filter  measurement  update  rate  is  only  0.1  Hz  (i.e.  every  ten  seconds). 
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Measurements  Based  On  Doppler  Radar  - 

Recall  that  the  fundamental  measurements  produced  by  the  Doppler  radar  system  are  the  body 
axis  velocity  components  Uq,  Vq,  Wo  (forward,  to  the  right  and  down  being  the  positive  senses  in  the 
body  axis  coordinate  frame).  It  is  then  necessary  to  process  velocity  differences  between  the  IRS  and  the 
Doppler  system  in  this  body  axis  frame.  In  order  to  do  this  properly,  the  IRS  velocity  components  must  be 
transformed  into  equivalent  body  axis  components  and,  as  well,  the  Doppler  velocity  components  must 
be  corrected  for  lever  arm  effects.  The  lever  arm  effects  are  due  to  the  fact  that  the  location  of  the  Doppler 
radar  antenna  is  not  coincident  with  the  location  of  the  IRS. 

Let  Cbc  be  the  transformation  matrix  that  converts  IRS  velocity  data  in  geographic  coordinates  (i.e. 
vn,  vB,  v2)  into  equivalent  components  in  body  axis  coordinates  (i.e.  U|,  V|,  W|).  This  transformation  matrix 
is  re-computed  continuously,  based  on  the  Euler  angle  (i.e.  attitude)  data  available  from  the  LTN-90-100 
IRS,  and  the  individual  matrix  element  calculations  are  shown  in  Appendix  C.  The  IRS  velocity 
components,  converted  into  body  axis  coordinates,  can  then  be  computed  as 

U|  =  Cij'Vn  +  C12VE  +  Ci3Vz 

V |  =  C2,i'vn  +  C2,2  vE  +  C-2,3  V2 

W,  =  C31  Vn  +  C32  Ve  +  C33  Vz  (7) 


In  order  to  look  at  lever  arm  effects,  define  fci  as  the  lever  arm  position  vector  from  the  Doppler 
radar  antenna  tfi  the  LTN-90  IRS.  This  position  vector  is  measured  in  body  axis  coordinates,  and  has 
components  /x,  ly  and  4.  Define  cob  as  the  angular  body  rate  vector  of  the  IRS,  as  expressed  in  body  axis 
coordinates.  It  is  directly  available  from  the  LTN-90-100  dataport  and  has  components  coBx,  (oBy,  coBz. 
Corrected  Doppler  velocity  components  (i.e.  Udc.  Vqc.  Wqc)  are  then  calculated  as  (see  Appendix  C  for 
details) 


Udc  = 

UD  - 

“Bz  /y 

+  03By/z 

VDC  * 

+ 

Q 

> 

“Bz 

-  (Ob  x/z 

Wdc  = 

WD  - 

U&ylx 

+  (OBx  /y 

(8) 

The  three  Doppler-based  measurement  components  available  at  discrete  times  tk,  Z3(k),  z 4(k)  and 
Z5(k),  can  then  be  written  as 

za(k)  =  U 1  (k)  -  UDC(k) 

Z4(k)  =  V 1  (k)  -  VDC(k) 

zs(k)  =  W,(k)  -  WDC(k)  (9) 
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The  24-element  rows  of  H  that  are  associated  with  these  three  measurements  (i.e.  rows  3, 4  and 
5)  have  the  following  forms: 

H3  =  (0  0  0  h3  4  h3  5  fl3  6  0  .  .  0  h3  20  0  0  Il3  23  ^>3,24) 

H4  =  (0  0  0  h44  h4  5  h4  6  0 ...  0  h4  21  0  h4  23  (1424) 

H5  =  (0  0  0  h54  hs  5  h5  6  0 . 0  hs  22  0  0)  (10) 


where  the  indicated  non-zero  elements  are  defined  in  Appendix  C. 

The  noise  variances  assumed  for  these  three  measurements,  r3,  r4  and  rs,  comprise  diagonal 
elements  3, 4  and  5  of  the  5  x  5  measurement  noise  covariance  matrix  R.  Nominal  values  chosen  for  them 
are  r3  =  (5.34  ft/s)2,  r4  =  (10.68  ft/s)2  and  rs  =  (5.34  ft/s)2,  based  on  an  analysis  of  typical  Doppler  radar  data 
from  the  Twin  Otter.  Doppler-based  measurements  are  available  at  a  2  Hz  rate,  much  faster  than  the  basic 
Kalman  filter  update  rate  of  0.1  Hz. 

Due  to  the  basic  nature  of  Bierman's  UDUT  factorization  algorithms,  the  individual  Loran  and 
Doppler-based  measurement  components  and  their  associated  statistics  (i.e.  z,  and  Hj,  r  j)  are  processed 
in  a  one-at-a-time  fashion.  Recall  that  the  software  algorithms  for  doing  this  form  of  Kalman  filter 
processing  have  been  described  in  subsection  2.1 . 


3.4  Measurement  Averaging  (Prefiltering} 

Kalman  filter  measurement  averaging  (or  prefiltering,  as  it  is  sometimes  called  -  see  Ref.  8)  is  a 
technique  used  to  take  full  advantage  of  discrete  measurement  data,  even  when  it  is  available  at  a  much 
higher  sample  rate  than  the  basic  update  rate  of  the  Kalman  filter  itself.  Rather  than  discarding  all  of  the 
measurement  information  'in  between'  Kalman  filter  time  update  points,  every  measurement  data  point  is 
used.  The  mathematical  details  of  this  concept  are  given  in  Appendix  D.  However,  a  short  description  of 
the  measurement  averaging  concept  is  given  below. 

The  use  of  the  measurement  averaging  concept  affects  the  way  in  which  z,  H,  v  and  R  are 
computed,  as  well  as  the  order  in  which  the  five  Kalman  filter  recursion  equations  (i.e.  Eqns.  3)  are 
invoked.  Let  AT  =  tk+i  -  tk  be  the  (constant)  Kalman  filter  update  interval,  and  assume  that  the  vector 
measurements  occur  at  N  equally  spaced  time  points  within  AT  (i.e.  z(tk  +  8tj)  where  5tj  =  i  AT/N;  i  = 
1,2,..,N).  Then  define  the  following  averaged  quantities: 

N 

zav  (k)  =  (1/N)  Iz(tk  +  5tj) 
i=1 

N 

Hav(k)  =  (1/N)  XlH(tk  +  8tj)  <*tk.tk  + S*i)] 
i=1 

N 

Vav(k)  =  (1/N)  Xv(tk  +  Sti) 
i=1 


N 

Rav(k)  =  (1/N2)  IR(tk+aj)  =  (1/N)  R 
i=1 
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(11) 


where  zav (k)  is  the  averaged  measurement,  Hav(k)  is  the  averaged  observation  matrix,  0(tk,  tk  +  8tj)  is  the 
so-called  intermediate  transition  matrix,  vav(k)  is  the  equivalent  averaged  measurement  noise  and  Rav(k)  is 
the  averaged  measurement  noise  covariance  matrix.  All  of  these  averaged  quantities  are  referenced  in 
time  to  tk,  and  Rav(k)  can  be  simplified  as  shown  above  when  it  is  assumed  to  be  a  constant  matrix.  With 
the  foregoing  definitions,  the  equivalent  averaged  observation  matrix  is 


zav(k)  =  H*(k)  8x(k)  +  vav(k)  (12) 


Kalman  Filter  Sequence  For  Prefllterfna  - 

The  Kalman  filter  equation  sequence,  as  shown  in  Eqns.  3,  must  be  modified  when  prefiltering 
takes  place.  The  filter  sequence  with  prefiltering  proceeds  as  follows: 

I)  At  current  time  tk+r  (i.e.  the  time  just  prior  to  processing  the  latest  averaged  measurements)  the 
Kalman  filter  acquires  an  averaged  measurement  vector  zav(k)  from  the  prefiltering  process,  with  zav(k) 
computed  by  averaging  the  measurements  sampled  over  the  update  interval  tk  -->  tk+i . 

ii)  The  Kalman  filter  processes  each  element  of  zav  sequentially  (i.e.  as  a  scalar  measurement),  using 
the  UDUPDC  subroutine,  in  order  to  compute  measurement  updates,  8x  and  P,  referenced  to  time  tk+. 

iii)  The  Kalman  filter  estimates,  8x(k+)  and  P(k+),  are  then  time-extrapolated  to  the  current  time  tk+i 
using  subroutine  MWGSUD.  The  time  update  of  8x,  as  computed  by  MWGSUD,  is  equivalent  to  the 
equation 


8x(k+T)  =  <t(k,k+1)  8x(k+) 


(13) 


Estimates  of  the  error  state  vector  8x(k+T)  are  then  in  time  synchronization  with  the  latest  values  of  the 
navigation  variables  available  from  any  of  the  navaids. 

iv)  In  the  case  of  the  assumed  feedforward  Kalman  filter  error  state  configuration,  8x(k+1)  contains 
the  cumulative  error  build-up  in  the  IRS  parameters  of  interest  at  tk+i  •  These  IRS  parameters  can  be 
corrected  at  the  basic  filter  rate  of  0.1  Hz,  or  even  at  the  higher  IRS  rate  of  16  Hz  (i.e.  by  applying  the  last 
best  estimates  of  8x  over  the  entire  10  second  filter  update  interval).  If  the  10  second  'resets'  in  the 
corrected  16  Hz  IRS  digital  parameters  appear  too  prominently,  there  is  even  the  possibility  of  low-pass 
filtering  (i.e.  digitally)  a  1 6  Hz  version  of  the  correction  signal  in  order  to  smooth  out  the  high  frequency 
component  of  the  resets. 

3.5  Euler  Angle  Error  Correction 

The  Kalman  filter  error  state  estimates  for  IRS  position  and  velocity  errors  are  computed  in  an 
appropriate  geographic  coordinate  frame  for  direct  correction  of  the  corresponding  IRS  parameters. 
However,  in  the  case  of  the  IRS  attitude  errors  (i.e.  Euler  angle  errors)  a  corrector  algorithm  must  be  used 
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in  order  to  transform  Kalman  filter  estimates  of  tilt  errors  into  equivalent  estimates  of  Euler  angle  errors. 
The  fundamental  tilt  errors  being  estimated  by  the  Kalman  filter,  e  =  [cm  eg  ej1,  are  the  attitude  errors  of 
the  so-called  computed  geographic  coordinate  frame  (i.e.  local  level,  N-E-z)  relative  to  the  true  geographic 
frame.  A  corrector  algorithm  converts  tilt  errors  into  equivalent  errors  in  the  direction  cosine  matrix  (DCM) 
that  describes  the  transformation  from  body  axes  to  geographic  axes,  denoted  by  CbG.  Once  the  DCM 
has  been  corrected,  it  is  quite  straightforward  to  compute  corrected  values  of  the  associated  Euler  angles 
(i.e.  0,  <|>,  y). 

Consider  e  to  be  a  small  angle  error  vector  representing  the  small  rotation  of  the  computed 
geographic  frame  about  the  true  geographic  (N,  E,  z)  axes.  Define  ex  as  the  matrix-equivalent  of  the  so- 
called  skew  symmetric  cross  product  operator  (Ref.  13)  that  can  be  associated  with  e.  Then  ex  has  the 
following  form: 


r 

0 

•Ez 

eg 

1 

ex 

Ez 

0 

-EN 

L 

-eg 

EN 

0 

J 

(14) 

Let  CbG*  be  the  error-corrupted  version  of  the  DCM  which  results  from  errors  in  the  Euler  angles 
being  measured  by  the  IRS.  in  Appendix  A  it  is  shown  how  CbG*  is  computed  using  the  IRS-supplied 
Euler  angles  0,  <}>  and  y.  Let  CbcG  be  a  corrected  version  of  this  DCM,  based  on  knowledge  of  e.  Then 
the  error  in  CbG*.  5CbG,  would  be  defined  as 


5CbG  =  CbG'  -  CbG 


(15) 


where  CbG  is  the  true  DCM.  It  can  be  shown  (Ref.  13)  that,  to  first  order, 


8CbG  =  -(ex)  QjG 


(16) 


Then  a  corrected  version  of  CbG\  CbcG,  can  be  calculated  as 


CbcG  =  0,0*-  6CbG 


=  CfeG*  +  (ex)CbG 
s  Cb®  +  (ex)  CbG* 

s  (I  +ex)  Cb^  (17) 


to  first  order. 


17 


Once  CbG*  has  been  corrected,  as  above,  corrected  values  for  the  Euler  angles  can  be  computed 
as 


0c  =  arc  tan  [C^icAl -Ca.ic2)172] 

<te  =  arc  tan  [  C3>2c/C3,3c  ] 

Vc  =  arctan[C2,ic/Ci,iC  ]  (18) 


From  Eqns.  18,  it  should  be  noted  that  only  five  of  the  nine  DCM  elements  in  CbG  need  to  be  corrected  in 
order  to  correct  the  Euler  angles,  namely  Ci  i ,  C2,i.  C3,i ,  C3  2  and  C33. 


Computation  of  Euler  Angle  Error  Bounds  • 

As  well  as  the  calculation  of  0c,  <t>c.  Vc.  using  Eqns.  18  above,  it  would  be  useful  to  have  an 
equivalent  1 -sigma  uncertainty  bound  for  each  of  the  Euler  angle  calculations,  just  as  we  do  for  the  other 
Kalman  filter  error  state  estimates.  This  is  not  straightforward,  because  errors  in  the  Euler  angles  do  nol 
have  a  direct,  one-to-one  correspondence  with  the  tilt  error  states  (for  which  the  Kalman  filter  does 
compute  an  updated  covariance).  However,  an  ad  hoc  procedure  that  should  give  a  representative 
uncertainty  bound  for  each  of  the  Euler  angles  would  be  as  follows: 

1 )  For  each  of  the  three  tilt  error  bounds,  set  the  corresponding  tilt  error  component  to  that  value 
(with  the  other  two  components  set  to  zero)  and  calculate  a  value  of  CbcG  using  Eqn.  17.  Repeat  the 
procedure  using  the  negative  of  the  bound  in  each  case.  The  result  will  be  six  different  values  of  CbcG- 

2)  For  each  of  the  six  CbcG's,  compute  'corrected'  values  for  the  three  Euler  angles  using  Eqn.  18. 
There  will  be  six  values  for  each  of  the  three  Euler  angles. 

3)  For  each  'corrected'  value,  compute  an  Euler  angle  error  by  differencing  the  'corrected'  value  with 
the  corresponding  IRS-based  Euler  angle.  The  result  will  be  six  different  error  values  for  each  of  the  three 
Euler  angles. 

4)  For  each  of  the  three  Euler  angles,  do  an  RSS  (i.e.  square  root  of  the  sum  of  the  squares)  of  the 
six  values  in  order  to  arrive  at  a  representative  value  for  the  bound. 


Computation  of  True  Tilt  Errors  - 

When  using  simulated  flight  test  data  (thus  having  access  to  the  true  aircraft  trajectory  and  true 
navigation  system  errors  being  simulated),  it  would  be  beneficial  to  be  able  to  derive  the  true  IRS  tilt  errors 
from  knowledge  of  the  Euler  angle  errors.  Equations  15  and  16,  when  combined,  yield 

(ex)  CbG  -  Q,G  -  CbG* 
or 


(ex)  CbG  (CbG)-i  =  I  -  OF  (Cb®)1 
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or 


ex  =  I  -  Cb<3*  (CbG)T,  since  (CbG)-i  =  (CbG)T  (19) 

Note  that  CbG*  is  the  DCM  that  would  be  calculated  from  error-corrupted  simulated  IRS  attitude  data.  On 
the  other  hand,  CbG  would  be  calculated  using  the  reference  trajectory  attitude  (i.e.  the  true  attitude). 
Once  the  skew  symmetric  matrix  (ex)  is  calculated  from  Eqn.19,  the  matrix  elements  that  correspond  to  the 
tilt  errors  eg  and  ez  can  be  identified  and  compared  with  Kalman  filter  estimates  of  the  same  tilt  error 
parameters. 


4.0  SOFTWARE  FOR  NAVIGATION  SIMULATION  AND  KALMAN  FILTERING 

In  order  to  analyze  and  refine  proposed  Kalman  filter  integrated  navigation  designs  properly, 
accurate  simulation  studies  must  first  be  conducted.  The  heart  of  these  simulation  studies  is  a  set  of 
algorithms  for  generating  'real  world'  data  for  the  various  navigation  sensors,  based  upon  a  specified 
aircraft  flight  profile  and  known  error  statistics  for  each  sensor.  A  comprehensive  set  of  navigation  sensor 
simulation  software  has  been  developed  at  the  Flight  Research  Laboratory  to  model,  quite  accurately, 
realistic  aircraft  flight  trajectory  dynamics,  and  the  error  dynamics  of  such  airborne  navigation  sensors  as 
medium  accuracy  INS’s,  Doppler  velocity  sensing  (DVS)  systems  and  various  radio  position  fixing/ranging 
systems.  This  section  describes  a  set  of  simulation  software  modules  that  have  been  written  for  execution 
under  the  IBM  VM  Operating  System  resident  on  NRC’s  IBM  3090  mainframe  computer.  Results  will  be 
given  to  show,  for  a  typical  simulated  flight  trajectory,  how  realistic  navigation  sensor  data  are  created.  A 
general  description  will  be  given  of  the  comprehensive  Kalman  filtering  software  package  that  has  been 
developed  to  run  on  either  simulated  or  real  navigation  data  from  the  Twin  Otter's  suite  of  navigation 
sensors.  A  series  of  plots  derived  from  the  Kalman  filter  running  on  simulated  data  will  be  shown  to 
document  the  performance  and  accuracy  of  the  filter  design.  Finally,  the  section  will  conclude  with  a  list  of 
the  general  findings  from  the  extensive  series  of  simulation  studies  that  were  conducted. 

4.1  Software  for  Generating  Reference  Trajectory  Profile  Data 

a)  PROFST 

A  VM  Executive  routine,  PROFST  EXEC,  is  used  to  run  a  trajectory  profile  generation 
FORTRAN  program  with  the  FORTRAN  name  PROFST.  PROFST  creates  3-dimensional  inertial  flight 
trajectory  data  (in  metric  units)  at  specified  points  in  time  using  a  user-supplied  description  of  the  desired 
profile.  The  resulting  trajectory  profile  dataset  is  then  used  as  the  reference  for  the  true  aircraft  trajectory 
when  analyzing  Kalman  filter  accuracy  and  it  is  also  the  starting  point  for  creating  realistic,  error-corrupted 
navigation  data  from  the  various  navaids  being  simulated.  This  program  calls  a  series  of  general-purpose 
profile  generation  subroutines  (i.e.  PROGEN,  PROGNI,  SEVAL,  SPLINE),  which  are  based  on  the 
use  of  cubic  splines.  PROFST  takes  the  inertial  velocity  data  supplied  by  the  PROGEN  subroutine  and 
integrates  it,  using  Simpson’s  rule  of  integration,  to  generate  3-dimensional  flight  trajectory  position  data. 
PROFST  also  creates  appropriate  output  files  containing  all  of  the  flight  trajectory  parameters  of  interest. 
When  running  the  program,  the  user  specifies  the  desired  sample  rate  (in  Hz)  of  the  trajectory  data  to  be 
created  and  the  total  time  length  (in  seconds)  of  the  trajectory.  The  input  and  output  files  involved  in  the 
execution  of  PROFST  have  the  following  characteristics: 
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Input  File:  A  dataset  with  the  generic  name  PROGD  DAT  An  (where  n  is  a  dataset  number 

to  denote  the  particular  type  of  trajectory  file  being  created)  is  required.  PROGD 
DATAn  specifies  the  starting  3-dimensional  position  (Lat,  Long,  altitude)  for  the 
trajectory  in  its  first  record,  and  then  consists  of  records  of  desired  aircraft  dynamic  state 
(i.e.  attitude  and  speed)  at  selected  points  along  the  trajectory,  to  a  maximum  of  100 
trajectory  points.  During  portions  of  the  trajectory  where  the  aircraft  state  changes  rapidly 
the  trajectory  points  are  more  closely  spaced  than  during  'benign'  sections  of  the 
trajectory.  This  approach  ensures  that  a  realistic  trajectory  is  'captured'  by  the  cubic  spline 
fitting  process  during  the  high  dynamics  portions  of  the  simulated  flight  trajectory. 

Output  Files:  PROFST  takes  the  sparse  set  of  trajectory  points  defined  in  the  input  dataset 
PROGD  DATAn  and  interpolates  between  the  points  using  cubic  spline  techniques. 
Two  output  files  are  normally  created  from  a  PROFST  run,  as  follows: 

SPUNFMT  OUT33  -  contains  20  flight  parameters  at  a  64  Hz  data  rate  (unformatted, 
single  precision),  with  64  Hz  chosen  as  the  minimum  rate  possible  in  order  to 
generate  strapdown  IMU  data  with  sufficient  accuracy.  After  the  idealized  strapdown 
IMU  data  is  generated,  this  rather  large  file  is  no  longer  needed  and  is  erased  to 
save  disk  storage  space. 

SPUNFMT  TWOHZ  -  contains  the  same  20  flight  parameters  as  SPUNFMT  OUT33 
(unformatted,  single  precision),  but  at  the  lower  rate  of  2  Hz.  This  file  is  a  lot  smaller 
than  SPUNFMT  OUT33  and  is  retained  for  use  as  the  reference  trajectory  data  set 
for  subsequent  Kalman  filter  analysis  purposes. 

Both  output  files,  SPUNFMT  OUT33  and  SPUNFMT  TWOHZ,  contain  the  same  set  of  20 
aircraft  inertial  flight  trajectory  parameters,  but  at  the  differing  data  rates.  For  each  dataset  the  20 
parameters  are  written  in  the  following  order: 

-  Trajectory  point  number  (i.e.  record  counter) 

-  Time  from  start  of  flight 

-  3  components  of  aircraft  velocity  (north,  east,  down  in  m/s) 

-  3  components  of  aircraft  velocity  rate  (north,  east,  down  in  m/s2) 

-  3  components  of  aircraft  attitude  (roll,  pitch,  heading  in  rad) 

-  3  components  of  aircraft  attitude  rate  (roll  rate,  pitch  rate,  yaw  rate  in  rad/s) 

-  3  components  of  position  change  from  start  location  (ALat,  ALong,  Aalt  in  deg,  deg,  m) 

-  3  components  of  absolute  aircraft  position  (Lat,  Long,  alt  in  deg,  deg,  m) 

b)  PROGPLT 

PROGPLT  EXEC  is  a  menu-driven  Executive  routine  to  plot  various  subsets  of  trajectory  profile 
parameters  that  have  been  created  by  a  PROFST  run  that  produces  a  SPUNFMT  OUT33  data  set. 
PROGPLT  calls  a  general-purpose,  FORTRAN-based  time  series  plotting  program,  PLOTTER  (with 
associated  subroutines  PLOTT  and  BGPLT),  which  has  been  written  at  the  Flight  Research  Laboratory 
and  utilizes  the  DISSPLA  package  of  plotting  subroutines.  There  are  five  different  subsets  of  flight 
trajectory  parameters  that  can  be  plotted,  namely  i)  the  3  inertial  velocity  components  (note  that  metric 
units  are  used),  ii)  the  3  aircraft  velocity  rate  components,  iii)  the  3  aircraft  attitude  components  (i.e.  Euler 
angles),  iv)  the  3  attitude  rate  components  and  v)  the  3  change-in-position  components.  Figures  4  to  8 
show  a  series  of  plots  generated  by  PROGPLT  for  a  1  hour  simulated  flight  profile  corresponding  to  a 
square  pattern  with  take-off  and  landing  included  as  well. 
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4.2  Software  for  Generating  Realistic  Navigation  Data 

A  series  of  programs  has  been  written  to  take  the  ideal  flight  trajectory  data,  model  realistic  levels 
of  time-varying  errors  in  each  of  the  navaids  to  be  simulated,  and  generate  data  files  of  simulated 
navigation  data  in  formats  similar  to  real  navigation  data  that  would  be  acquired  onboard  the  aircraft.  At  the 
moment,  the  software  is  configured  to  handle  generic  versions  of  a  strapdown  INS,  Doppler  velocity 
sensor  and  radio  position  fixing  sensor.  For  this  report,  the  software  has  been  used  to  accurately  simulate 
an  LTN-90-100  IRS,  a  Decca  Type  72  Doppler  radar  and  an  ARNAV  R-40  Loran-C  receiver. 

a)  SDIMU 

By  far  the  most  complex  navigation  simulation  software  to  develop  is  that  required  to  accurately 
simulate  an  IMU  and/or  INS.  SDIMU  EXEC  is  an  Executive  routine  that  is  used  to  run  SDIMU,  a 
FORTRAN  program  that  simulates  ideal  (i.e.  perfectly  accurate)  data  being  generated  from  a  strapdown 
IMU,  as  well  as  ideal  height  (i.e.  altimetry)  information.  The  raw  outputs  from  a  strapdown  IMU  are  the  body 
frame  velocity  and  attitude  increments  (i.e.  changes  in  velocity  and  attitude  over  a  very  small  time  interval) 
usually  sampled  at  a  fairly  high  data  rate.  For  simulation  purposes,  IMU  increments  are  calculated  at  a  64  Hz 
rate  by  running  SDIMU  using  the  inertial  rate  information  available  from  the  SPUNFMT  OUT33 
reference  trajectory  data  file  as  input.  A  small  input  control  file,  SDIMU  CTRL,  specifies  the  strapdown 
IMU  sample  rate  and  the  time  length  desired  for  the  run.  The  standard  output  from  SDIMU  is  another  64 
Hz  single  precision,  unformatted  data  file,  SDIMU  DATA,  which  consists  of  8  parameters,  namely  i)  a 
record  counter,  ii)  3  components  of  incremental  change  in  body  axis  velocity  (m/s),  iii)  3  components  of 
incremental  change  in  body  axis  attitude  (rad)  and  iv)  true  altitude  (m).  As  previously  stated,  upon 
completion  of  a  SDIMU  run,  the  64  Hz  input  file,  SPUNFMT  OUT33,  which  is  not  required  for 
subsequent  analysis  purposes,  is  erased  to  save  disk  storage  space. 


b)  SDNAV1 

An  Executive  routine,  SDNAVG1  EXEC,  controls  the  execution  of  SDNAV1 ,  a  FORTRAN 
program  that  processes  the  strapdown  IMU  incremental  angles  and  velocities  in  order  to  emulate  a 
strapdown  INS  or  IRS  (such  as  the  LTN-90-100)  computing  inertial  position,  velocity  and  attitude.  In 
SDNAV1 ,  the  navigation  equations  are  solved  using  a  wander  azimuth  navigator  mechanization,  with  all 
updates  occurring  at  a  64  Hz  rate  (as  is  done  in  the  LTN-90-100).  However,  the  standard  output  rate  for 
the  simulated  INS  data  is  16  Hz  which  is  the  rate  at  which  LTN-90-100  data  are  sampled  onboard  the  Twin 
Otter.  Included  in  SDNAV1  is  code  for  modelling  realistic  errors  in  a  baro  altimeter  source  as  well  as  an 
emulation  of  the  third-order  baro  damping  loop  that  is  implemented  inside  the  LTN-90-100.  The 
SDNAVl  software  was  developed  at  Flight  Research  as  a  major  re-structuring  of  a  DREO  software 
package  called  NAVGTR  and  utilizes  DREO-developed  subroutines  QCINI,  CORIOL,  QCUPDT  and 
NAVCAL  for  performing  the  various  standard  strapdown  navigation  calculations.  There  are  2  input  files 
required  in  order  to  run  SDNAVl  -  a  control  file,  SDNAVl  CTRL,  and  the  data  file,  SDIMU  DATA,  of 
ideal  strapdown  IMU  increments  sampled  at  64  Hz.  SDNAVl  CTRL  contains  the  following  information:  i) 
input/output  sample  rates,  ii)  time  length  of  the  run,  iii)  initial  conditions  on  aircraft  position,  velocity  and 
attitude  and  iv)  statistical  information  on  gyro,  accelerometer  and  altimeter  errors.  Initial  conditions  can  be 
specified  in  this  input  file  to  simulate  typical  INS  start-up  errors  arising  from  the  alignment  process.  The 
statistical  information  required  to  model  gyro  and  accelerometer  errors  properly  is  the  kind  shown  in  Table 
3.2  of  Section  3  for  the  gyros  and  accelerometers  in  the  LTN-90-100.  A  combination  of  bias,  ZMWG  noise 
and  first-order  Markov  error  modelling  is  used  in  SDNAVl  to  produce  realistic  time  varying  errors  in  the 
gyros,  accelerometers  and  altimeter.  The  standard  output  file  created  by  running  SDNAVl  is 
TWOTTER  INSDATA,  a  data  file  consisting  of  16  INS  parameters  sampled  at  16  Hz  in  the  same  order, 
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and  with  the  same  engineering  units,  as  real  LTN-90-100  IRS  data  that  would  be  collected  onboard  the 
Twin  Otter  (i.e.  as  shown  in  Table  3.1  of  Section  3). 


C)  SDERKF 

A  FORTRAN  program,  SDERKF,  is  used  to  compute  (for  plotting  purposes)  the  differences 
between  the  simulated  strapdown  INS  parameters,  as  created  by  SDNAV1  in  its  output  data  set 
TWOTTER  INSDATA,  and  the  original  inertial  reference  parameters  created  by  PROFST  in  its  output 
file  SPUNFMT  TWOHZ.  The  execution  of  SDERKF  is  normally  controlled  by  the  same  Executive 
routine  that  is  used  to  control  SDNAV1 ,  namely  SDNAVG1  EXEC  -  i.e.  a  run  of  SDNAV1  is 
immediately  followed  by  a  run  of  SDERKF  via  SDNAVG1  EXEC.  SDERKF  requires  3  input  files  in 
order  to  run:  SDNAV1  CTRL,  TWOTTER  INSDATA  and  SPUNFMT  TWOHZ.  It  produces  one 
main  output  file,  SPUNFMT  ERR,  which  consists  of  a  record  counter  plus  the  9  true  inertial  error 
quantities  corresponding  to  the  simulated  IRS  system  (3  position  component  errors,  3  velocity 
component  errors  and  3  attitude  component  errors)  all  sampled  at  a  0.1  Hz  rate.  The  engineering  units 
chosen  for  the  inertial  errors  are  the  same  ones  being  used  for  Kalman  filter  output  display  purposes  (see 
subsection  4.4). 


d)  ERRPLOT 

ERRPLOT  EXEC  is  a  menu-driven  Executive  routine  for  controlling  the  plotting  of  various 
subsets  of  simulated  IRS  errors  using  the  Flight  Research  general-purpose  time  series  plotting  package, 
PLOTTER.  The  user  can  chose  any  of  the  following  three  subsets  of  true  IRS  errors  to  be  plotted:  i)  the 
3  IRS  position  error  components,  ii)  the  3  IRS  velocity  error  components  or  iii)  the  3  IRS  attitude  error 
components.  Figures  9  to  1 1  show  typical  plotted  outputs  for  each  of  these  subsets.  In  this  particular 
case,  the  flight  trajectory  consists  of  a  1  hour  square  pattern,  including  takeoff  and  landing  phases. 
ERRPLOT  has  been  used  to  confirm  the  accuracy  of  the  IRS  simulation  software  by  comparing  plotted 
outputs  with  those  for  known  IRS  behaviour  (Refs.  14,  15)  under  the  same  set  of  initial  conditions  and 
sensor  errors. 


e)  DOPLOR 

DOPLOR  EXEC  is  used  to  control  the  execution  of  DOPLOR,  a  FORTRAN-based  program 
that  creates  realistic  simulated  navigation  data,  based  on  a  desired  reference  trajectory,  that  corresponds 
to  data  from  a  Decca  Type  72  Doppler  radar  and  an  ARNAV  R-40  airborne  Loran-C  receiver,  including  data 
rates  and  engineering  units  identical  to  those  for  real  data  from  the  Twin  Otter  data  acquisition  system. 
DOPLOR  requires  2  input  files  in  order  to  run  -  a  control  file,  DOPLOR  CTRL,  and  the  2  Hz  reference 
trajectory  data  file,  SPUNFMT  TWOHZ,  that  has  been  created  by  PROFST.  DOPLOR  CTRL 
contains  the  following  information  required  for  the  proper  simulation  of  Doppler  and  Loran-C  data:  i)  time 
length  of  the  run  and  time  interval  for  overwater  flight,  ii)  Doppler-to-INS  lever  arms,  iii)  Loran-C  latitude  and 
longitude  error  statistics  and  iv)  Doppler  error  statistics.  Loran-C  error  statistics  consist  of  ZMWG  noise  and 
first-order  Markov  specifications,  while  Doppler  error  statistics  consist  of  first-order  Markov  processes  for 
each  of  U  channel  scale  factor,  V  and  W  channel  boresights  and  ZMWG  noise  in  each  of  the  U,  V  and  W 
channels.  As  well,  if  overwater  flight  is  involved,  then  first-order  Markov  processes  are  defined  for  the  N 
and  E  sea  bias  components.  Typical  values  for  these  various  statistics  can  be  found  in  Tables  3.5  and  3.6 
and  sub-subsection  3.3.3  of  Section  3.  Two  output  files  are  created  from  a  DOPLOR  run  -  i)  TWOTTER 
LORDATA  which  contains  realistic,  error-corrupted  Loran-C  Lat  and  Long  data  and  the  true  Loran-C 
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Markov  error  processes,  all  at  the  standard  Loran-C  sample  rate  of  1  Hz  and  ii)  TWOTTER  DOPDATA, 
containing  realistic,  error-corrupted  Doppler  U,  V,  W  channel  data  and  the  true  Doppler  Markov  error 
processes,  all  at  the  standard  Decca  Doppler  sample  rate  of  2  Hz.  The  true  Markov  error  processes  for  the 
simulated  Loran-C  and  Doppler  data  become  useful  for  assessing  the  accuracy  of  a  Kalman  filter  algorithm 
when  it  is  trying  to  estimate  these  particular  Markov  error  states. 


4.3  Kalman  Filter  Software  Descriptions 

A  comprehensive  set  of  Kalman  filter  software  has  been  developed  at  Flight  Research,  based  on 
the  Kalman  filter  design  described  in  Section  3,  for  use  with  Twin  Otter  navigation  data  in  either  real  or 
simulated  form.  The  purpose  of  this  subsection  is  to  outline  the  algorithmic  and  hierarchical  structure  of 
the  Kalman  filter  software  down  to  the  level  of  the  individual  subroutines  that  are  called  and  the 
input/output  datasets  that  are  normally  required. 


a)  SDKFNAV 

SDKFNAV  EXEC  is  a  menu-driven  Executive  routine  for  controlling  the  execution  of  four 
slightly  different  versions  of  the  Kalman  filter  software.  These  four  different  versions,  called  by 
SDKFNAV  EXEC  under  user  control,  have  the  following  main  FORTRAN  program  names:  i)  SDKFN1  - 
for  running  the  Kalman  filter  on  real  data  with  the  measurement  averaging  option  (see  subsection  3.4  for 
details),  ii)  SDKFN2  -  for  running  the  Kalman  filter  on  real  data,  but  with  no  measurement  averaging,  iii) 
SDKFN3  -  for  running  the  Kalman  filter  on  simulated  data  with  the  measurement  averaging  option  and  iv) 
SDKFN4  -  for  running  the  Kalman  filter  on  simulated  data  with  no  measurement  averaging.  Within 
SDKFNAV  EXEC,  one  of  two  other  EXEC'S  is  invoked  depending  on  whether  real  or  simulated  data  is 
to  be  used  -  i)  SDKFFDEF  EXEC  defines  all  of  the  input  and  output  files  that  are  required  for  the 
simulated  data  case  and  ii)  SDKFRFDF  EXEC  defines  all  of  the  input  and  output  files  for  the  real  data 
case.  The  following  are  very  brief  descriptions  of  the  various  input  and  output  files  that  are  required  for  the 
various  versions  of  the  Kalman  filter  as  well  as  descriptions  of  the  complete  set  of  FORTRAN  modules  that 
comprise  the  Kalman  filter  software  package: 

Incut  Files: 

SDKFNAV  CTRL  -  passes,  to  any  version  of  SDKFNn,  all  of  the  control  information  that  is 

required  to  initialize  the  filter,  including:  i)  time  length  of  the  filter  run  and  time 
interval  for  overwater  flight,  ii)  nominal  starting  latitude  and  Doppler-to-INS 
lever  arms,  iii)  IRS  system  error  state  statistics  (i.e.  Table  3.4),  iv)  sensor 
Markov  error  statistics  (i.e.  Table  3.5),  v)  IRS  system  plant  noise  statistics  (i.e. 
Table  3.6)  and  vi)  measurement  noise  statistics  (given  in  sub-subsection 
3.3.3).  As  well,  for  the  real  data  case  there  may  be  further  information 
supplied  at  the  end  of  this  file  concerning  Visual  On-Top  (VOT)  position 
reference  points  along  the  flight  path. 

-  the  1  Hz  version  of  airborne  Loran-C  data,  having  the  identical  format  for 
both  the  real  and  simulated  data  cases,  except  that  in  the  simulated  data 
case  the  true  Loran-C  Markov  error  processes  are  added  on  to  the  end  of 
each  record  for  Kalman  filter  analysis  purposes. 


TWOTTER  LORDATA 
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TWOTTER  DOPDATA  - 

TWOTTER  INSDATA  - 

SPUNFMT  ERR  - 

Output  Files: 

TWOTTER  MEASURE  - 

TWOTTER  SYSER1  - 

TWOTTER  SYSER2  - 

TWOTTER  MARKER1  - 

TWOTTER  MARKER2  - 

TWOTTER  INERRS  - 

TWOTTER  LRERRS  - 

TWOTTER  KFERRS  - 


the  2  Hz  version  of  the  Decca  Doppler  radar  data,  having  identical  formats  for 
both  the  real  and  simulated  data  cases,  except  that  in  the  simulated  data 
case  the  true  Doppler  Markov  error  processes  are  added  on  to  the  end  of 
each  record  for  Kalman  filter  analysis  purposes. 

the  16  Hz  version  of  the  Litton  LTN-90-100  IRS  data,  as  defined  in  Table 
3.1 ,  having  identical  formats  for  both  the  real  and  simulated  data  cases. 

the  0.1  Hz  version  of  true  IRS  system  errors,  as  produced  by  program 
SDERKF,  and  used  by  the  Kalman  filter  program,  SDKFNn  (n  =  3  or  4),  to 
evaluate  its  state  estimation  accuracy.  This  input  file  is  only  used  for  the 
simulated  data  case. 


a  0.1  Hz  data  file  consisting  of  the  five  scalar  measurements  being 
processed  by  the  Kalman  filter,  namely  i)  2  measurements  based  on  Loran-C 
and  ii)  3  measurements  based  on  Doppler  radar  (see  sub-subsection  3.3.3 
for  details).  It  is  instructive  to  plot  this  measurement  data  to  confirm  that  it 
looks  reasonable  for  processing  by  the  Kalman  filter, 
a  0.1  Hz  unformatted  data  file  of  filter-estimated  IRS  system  errors  and 
associated  +/-  1 -sigma  bounds,  for  plotting  purposes.  This  particular  file 
contains  information  on  the  first  6  IRS  system  error  states,  as  shown  in  Table 

3.4. 

a  0.1  Hz  unformatted  data  file,  similar  to  TWOTTER  SYSER1,  but 
containing  information  on  the  last  4  IRS  system  error  states  shown  in  Table 

3.4,  for  plotting  purposes. 

a  0.1  Hz  unformatted  data  file  of  filter-estimated  sensor  Markov  error  states 
and  associated  +/- 1 -sigma  bounds,  for  plotting  purposes.  This  particular  file 
contains  information  on  the  first  7  Markov  error  states,  as  shown  in  Table  3.5. 

a  0.1  Hz  unformatted  data  file,  similar  to  TWOTTER  MARKER1 , 
containing  information  on  the  last  7  sensor  Markov  error  states  shown  in 
Table  3.5  for  plotting  purposes. 

a  small  file  containing  IRS  horizontal  position  error  information  at  each  of  the 
Visual  On-Top  (VOT)  points  overflown  during  an  actual  Twin  Otter  navigation 
flight.  This  output  file  is  only  created  for  the  real  data  case  and  then  used  to 
plot  IRS  position  errors. 

a  small  file  containing  Loran-C  horizontal  position  error  information  at  each  of 
the  VOT  points  overflown  during  an  actual  navigation  flight.  This  file  is  only 
created  for  the  real  data  case  and  used  for  plotting  actual  Loran-C  position 
errors  during  a  flight. 

a  small  file  containing  filter-corrected  IRS  horizontal  position  error  information 
at  each  of  the  VOT  points  overflown  during  an  actual  flight.  This  file  is  only 
created  for  the  real  data  case  and  used  for  determining  the  true  position 


accuracy  of  the  Kalman  filter. 


TWOTTER  KFBNDS  -  a  0.1  Hz  data  file  of  filter-generated  +/- 1  -sigma  bounds  for  the  IRS  horizontal 

position  error  states.  This  information,  only  created  in  the  real  data  case,  is 
plotted  along  with  that  of  TWOTTER  KFERRS  in  order  to  evaluate  Kalman 
filter  performance. 


SDKFNn  -  this  is  the  main  Kalman  filter  program,  from  which  the  various  subroutines  associated  with 
specific  filtering  tasks  are  called.  As  already  stated,  there  are  actually  4  slightly  different 
versions  of  the  Kalman  filter  software  (in  order  to  handle  the  options  of  real  or  simulated 
data  and  prefiltering  or  no  prefiltering),  with  most  of  the  program  coding  differences  taking 
place  in  the  4  versions  of  main  program  SDKFNn  (n  =  1 ,. . .  ,4).  The  main  program  sets 
up  the  specified  initial  conditions  for  the  filter  run,  controls  the  0.1  Hz  rate  at  which 
measurements  are  defined  and  the  various  Kalman  filter  updates  take  place  and  writes  out 
the  pertinent  data  to  the  various  output  files  that  are  created. 

SDPHGM  -  a  subroutine,  called  by  SDKFNn  at  a  2  Hz  rate,  that  calculates  the  time  varying  portion  of 
the  so-called  intermediate  <6  transition  matrix  (see  Appendix  D)  and  intermediate  G  matrix 
using  trapezoidal  integration  on  the  16  Hz  inertial  data.  A  slightly  different  version  of  this 
subroutine,  SDPHG1,  is  used  when  measurement  averaging  (i.e.  prefiltering)  is 
involved. 

SDFGCL  -  a  subroutine,  called  by  SDPHGM,  that  computes  the  individual  matrix  elements  of  F  and 
G  according  to  the  equations  given  in  Appendix  B.  As  well,  the  latest  elements  of  the 
DCM  (see  Appendix  A)  are  computed.  In  the  case  of  prefiltering,  a  slightly  modified 
version  of  this  subroutine,  SDFGC1,  is  used. 

SDDRCL  -  a  subroutine,  called  by  SDKFNn  at  a  2  Hz  rate,  that  calculates  the  IRS/Doppler  velocity 
difference  measurements,  including  lever  arm  corrections  (see  Appendix  C  for  details).  A 
slightly  modified  version  of  this  subroutine,  SDDCLC,  is  used  for  the  simulated  data 
case. 

SDHCAL  -  a  subroutine,  called  by  SDKFNn  at  a  0.1  Hz  rate,  that  calculates  the  non-zero  elements 
of  the  Hj  vectors  (i  =  3,  4,  5)  required  for  Doppler-based  measurement  processing.  A 
slightly  different  version  of  this  subroutine,  SDHCLC,  is  used  for  the  prefilter  case. 

FACTOR  -  a  general-purpose  Kalman  filter  subroutine,  called  by  SDKFNn  and  described  in 
subsection  2.2,  for  the  efficient  factorization  of  the  error  state  covariance  matrix  P  into 
Bierman's  UDUT  form. 

UNFACT  -  a  general-purpose  Kalman  filter  subroutine,  called  by  SDKFNn  and  described  in 
subsection  2.2,  for  reconstructing  a  P  matrix  from  its  UDUT  decomposition  form.  It  is 
used  to  derive  the  +/-1  -sigma  bounds  of  the  various  error  state  estimates,  at  each  Kalman 
filter  update,  for  display  purposes. 

MWGSUD  -  a  general-purpose  Kalman  filter  subroutine,  called  by  SDKFNn  and  described  in 
subsection  2.2,  for  performing  Bierman's  efficient  time  update  of  the  error  state  estimate 
x  and  its  factored  covariance  matrix  U. 
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UDUPDC  -  a  general-purpose  Kalman  filter  subroutine,  called  by  SDKFNn  and  described  in 
subsection  2.2,  for  performing  Bierman's  efficient  measurement  update  of  x  and  U  based 
on  individual  scalar  measurements. 


NMDLAT  -  a  small  subroutine  used  in  SDKFNn  for  accurately  converting  a  change  in  latitude,  as 
expressed  in  units  of  nm,  into  the  equivalent  change  in  units  of  geographical  degrees  by 
using  local  earth's  curvature  information. 

NMDLNG  -  a  small  subroutine  used  in  SDKFNn  for  accurately  converting  a  change  in  longitude,  as 
expressed  in  units  of  nm,  into  the  equivalent  change  in  units  of  geographical  degrees 
using  local  earth’s  curvature  information. 


CONLAT  -  a  small  subroutine  used  in  SDKFNn  for  accurately  converting  a  change  in  latitude,  as 
expressed  in  units  of  geographical  degrees,  into  the  equivalent  change  in  units  of  nm 
using  local  earth's  curvature  information. 


CONLNG  -  a  smail  subroutine  used  in  SDKFNn  for  accurately  converting  a  change  in  longitude,  as 
expressed  in  units  of  geographical  degrees,  into  the  equivalent  change  in  units  of  nm 
using  local  earth's  curvature  information. 


b)  MESPLT 

MESPLT  EXEC  is  a  menu-driven  Executive  routine  for  controlling  the  execution  of  the 
FORTRAN  plotting  program,  MESPLT,  that  plots  the  2  different  types  of  Kalman  filter  measurement  data 
being  processed,  namely  i)  the  Loran-C/IRS  position  difference  measurements  ( z ^  and  z2  in  sub¬ 
subsection  3.3.3)  and  ii)  the  Doppler/IRS  velocity  difference  measurements  (z3,  z4  and  z5  in  sub¬ 
subsection  3.3.3).  The  data  to  be  plotted  is  contained  in  file  TWOTTER  MEASURE  which  is  created 
by  SDKFNn. 

C)  VOTPLT 

VOTPLT  EXEC  is  a  menu-driven  Executive  routine  for  controlling  the  execution  of  the 
FORTRAN  program,  VOTPLT,  that  plots  any  or  all  of  the  following  sets  of  errors:  i)  Kalman  filter  horizontal 
position  errors  at  the  VOT  points,  ii)  IRS  horizontal  position  errors  at  the  VOT  points  and  iii)  Loran-C 
horizontal  position  errors  at  the  VOT  points.  VOTPLT  is  only  used  for  the  real  data  case  and  requires  the 
following  input  data  files  as  created  by  SDKFNn.  TWOTTER  KFERRS,  TWOTTER  KFBNDS, 
TWOTTER  INERRS  and  TWOTTER  LRERRS. 

d)  SDKFPLT 

SDKFPLT  EXEC  is  a  menu-driven  Executive  routine  for  controlling  the  plotting  of  various 
subsets  of  Kalman  filter-generated  error  state  estimates  and  associated  +/-  1  -sigma  bounds  using  the 
Flight  Research  general-purpose  time  series  plot  program,  PLOTTER.  The  following  subsets  of 
parameters  can  be  chosen  by  the  user  for  plotting: 

I)  filter  estimates  of  the  3  IRS  position  error  components  (or  actual  filter  errors  in  measuring  the  IRS 
position  error  states,  for  the  simulated  data  case)  plus  associated  +/- 1  -sigma  bounds, 
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ii)  fitter  estimates  of  the  3  IRS  velocity  error  components  (or  actual  filter  errors  in  measuring  the  IRS 
velocity  error  states,  for  the  simulated  data  case)  plus  associated  +/- 1 -sigma  bounds, 

Hi)  filter  estimates  of  the  3  IRS  tilt  error  components  (or  actual  filter  errors  in  measuring  the  tilt  errors, 
for  the  simulated  data  case)  and  associated  +/- 1 -sigma  bounds, 

tv)  filter  estimates  of  the  3  IRS  gyro  Markov  error  components  and  associated  +/- 1  -sigma  bounds, 

v)  filter  estimates  of  the  3  IRS  accelerometer  Markov  error  components  and  the  associated  +/- 1- 
sigma  bounds, 

vi)  filter  estimates  of  the  2  Loran-C  Markov  error  components  and  the  2  sea  bias  Markov  error 
components  (or  actual  filter  errors  in  estimating  these  components,  for  the  simulated  data  case)  plus 
associated  +/- 1 -sigma  bounds, 

vil)  filter  estimates  of  the  3  Doppler  Markov  error  components  (or  actual  filter  errors  in  estimating  these 
components,  for  the  simulated  data  case)  and  associated  +/- 1 -sigma  bounds, 

viii)  actual  Loran-C  Markov  error  processes  (for  the  simulated  data  case  only) 

lx)  actual  Doppler  Markov  error  processes  (for  the  simulated  data  case  only). 

The  following  datasets  are  required  in  order  to  invoke  all  options  of  SDKFPLT :  TWOTTER 
SYSER1,  TWOTTER  SYSER2,  TWOTTER  MARKER1 ,  TWOTTER  MARKER2,  TWOTTER 
LORDATA,  TWOTTER  DOPDATA. 

e)  INNPLT 

INNPLT  EXEC  is  a  menu-driven  Executive  routine  for  controlling  the  execution  of  the 
FORTRAN  plotting  program,  INNPLT,  that  plots  the  5  measurement  residual  time  series  (described  in 
subsection  4.4). 

4.4  Kalman  Filter  Results  Using  Simulated  Data 

Not  a  lot  of  emphasis  will  be  placed  on  these  simulation  results,  since  the  most  important  results 
are  certainly  those  demonstrating  Kalman  filter  performance  for  real  navigation  data  taken  under  typical 
operating  conditions.  However,  in  order  to  demonstrate  that  the  fundamental  Kalman  filter  design  is 
sound,  and  performs  according  to  the  theory,  some  simulated  data  results  will  be  shown.  This  will  also 
serve  as  a  demonstration  of  the  full  capabilities  of  the  navigation  simulation/Kalman  filtering  software 
package. 

Kalman  filter  results  will  he  shown  for  the  same  one  hour  simulated  flight  scenario  that  was  used  to 
show  example  output  plots  from  the  ERRPLOT  plotting  routine  (i.e.  a  flight  trajectory  consisting  of  a 
square  pattern  together  with  takeoff  and  landing  phases  -  see  Figs.  4  to  8).  Table  4.1  gives  the  details  of 
the  IRS  initialization  errors  and  sensor  statistical  errors  that  were  used  in  control  file  SDNAV1  CTRL 
when  running  SDNAV1  to  produce  the  simulated  IRS  navigation  data  having  the  time  varying  errors 
shown  in  Figs.  9  to  1 1 .  From  these  figures  it  can  be  seen  that  there  is  a  significant  build-up  in  IRS  latitude 
and  north  velocity  error  in  particular;  however,  these  and  the  other  IRS  error  levels  are  still  realistic  (both  in 
their  initial  values  and  variations  with  time),  based  on  the  specifications  given  by  Litton  for  the  LTN-90-100 
IRS. 
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Table  4.2  shows  the  statistical  error  specifications  used  for  simulating  Loran-C  and  Doppler 
navigation  data  for  the  same  flight  trajectory  as  defined  above  for  simulated  IRS  data.  Figures  12  to  13 
show  the  simulated  first-order  Markov  error  time  histories  that  define  the  time  varying  errors  in  each  of 
these  navaids.  The  plots  shown  in  these  figures  were  created  using  appropriate  options  within  the 
SDKFPLT  plotting  package  described  in  subsection  4.3.  The  fundamental  time  varying  nature  of  these 
Markov  simulations  has  been  compared  with  the  known  error  behaviour  of  the  actual  navigation  devices  to 
confirm  the  realism  of  the  simulation  software. 

Figure  14  is  actually  a  copy  of  the  SDKFNAV  CTRL  file  for  a  representative  one  hour  Kalman 
filter  run,  using  the  SDKFNAV  software  package,  in  order  to  demonstrate  the  theoretical  capabilities  of 
the  filter.  For  this  particular  run,  filter  specifications  of  the  error  statistics  for  the  various  navaids  (shown  in 
the  listing  of  SDKFNAV  CTRL  -  i.e.  Fig.  14)  were  matched  closely  to  the  actual  error  statistics  being 
used  to  generate  the  simulated  navigation  data.  When  there  is  a  perfect  match  between  actual  sensor 
error  statistics  and  the  Kalman  filter  specification  of  the  same  parameters,  the  filter  is  ideally  'tuned'  and 
should  exhibit  optimal  performance.  Figures  15  and  16  show  one  hour  duration  time  history  plots  of  the 
two  different  measurement  types  (a  total  of  5  scalar  measurements  available  at  each  filter  update  interval). 
Again,  these  measurement  process  time  histories,  based  on  simulated  navigation  data,  have  been 
compared  with  those  for  representative  real  navigation  data  to  confirm  their  accuracy.  The  high  frequency 
portion  of  each  scalar  measurement  time  history  should,  and  does,  conform  to  the  Kalman  filter 
specification  of  the  1 -sigma  level  of  ZMWG  measurement  noise  for  the  associated  measurement 
component. 


TABLE  4.1 

ERROR  SPECIFICATIONS  FOR  LTN-90-100  IRS  SIMULATION 


Initial  Position  Errors: 
Initial  Velocity  Errors: 
Initial  Attitude  Errors: 
Gyro  Biases: 

Gyro  Drifts: 


A  LAT  =  0.00167 
A  vn  =  0.33  Ft/S 
A  ROLL  =  0.012  c 
0.01  Deg/Hr 


ALONG  =  0.00833  ° 
A  Vg  =  0.33  Ft/S 
APCH  =  0.012° 
Gvro  Scale  Factors: 
Correlation  Times: 
0.001  Deg/Hr  (la) 
Accel  Scale  Factors: 
Correlation  Times: 
5.0  p  G  (la) 


A  ALT  =  3.3  Ft 
A  vz  =  0.33  Ft/S 
A  HDG  =  0.0167° 
5  x  10-6 
3600  Sec 

5x10-5 
7200  Sec 


0.003  Deg/Hr 

Gyro-Bandom  Noise  Components: 

Accel  Biases:  50  pG 

Accel  Drifts:  5.0  pG 

Accel  Random  Nuise  Components: 

Baro  Scale  Factor:  0.05  CorrelTime:  5000  Sec 


Random  Noise:  0.1  m(1a) 
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TABLE  4.2 

ERROR  SPECIFICATIONS  FOR  LORAN-C  AND  DOPPLER  SIMULATIONS 

LORAN-C  STATISTICS: 

Latitude  Offset:  0.1  Nm  Latitude  Markov:  0.2  Nm  Correl  Time:  7200  Sec 

Longitude  Offset:  0.1  Nm  Longitude  Markov:  0.2  Nm  Correl  Time:  7200  Sec 

Latitude  ZMWG  Noise:  0.05  Nm  (la)  Longitude  ZMWG  Noise:  0.05  Nm  (la) 


DOPPLER  STATISTICS: 


SFu  Offset:  0.01 
BORv  Offset:  0.5  Deg 
BORw  Offset:  0.5  Deg 
SBn  Offset:  2.0  Ft/S 
SBE  Offset:  2.0  Ft/S 
U  Ch  Noise:  5.34  Ft/S  (la) 


SFu  Markov:  0.01 
BORv  Markov:  0.65  Deg 
BORw  Markov:  0.65  Deg 
SBn  Markov:  4.5  Ft/S 
SBe  Markov:  4.5  Ft/S 

V  Ch  Noise:  5.34  Ft/S  (la) 


SFy  Correl  Time:  15,000  Sec 
BORy  Correl  Time:  15,000  Sec 
BORw  Correl  Time:  15,000  Sec 
SBn  Correl  Time:  900  Sec 
SBe  Correl  Time:  900  Sec 
WCh  Noise:  5.34  Ft/S  (la) 


Figures  17  to  23  show  time  history  plots  that  demonstrate  the  Kalman  filter’s  accuracy  in 
estimating  various  sets  of  the  error  states  being  modelled.  The  +/- 1 -sigma  bounds,  as  computed  by  the 
filter,  are  also  displayed  along  with  the  associated  filter  state  error  to  indicate  the  level  of  confidence  that 
the  filter  is  assigning  to  that  particular  state  estimate.  According  to  the  theory  of  Kalman  filtering,  when  the 
filter  is  properly  'tuned'  the  actual  state  estimate  errors  should  be  within  the  +/-  1 -sigma  bounds 
approximately  67  %  of  the  time.  As  can  be  seen  from  these  plots  (and  results  from  many  other  simulation 
runs  as  well),  this  is  generally  the  case  for  this  particular  filter  run  and  confirms  that  the  Kalman  filter 
software  is  fundamentally  sound.  It  should  be  noted  that  this  Kalman  filter  run  is  of  short  duration 
compared  to  the  length  of  a  typical  operational  flight  and  there  is  always  an  initial  filter  'transient'  time  period 
that  shows  up  in  the  filter  state  estimates  while  the  filter  is  settling  down  to  a  steady  state  condition.  Plots 
of  the  estimates  of  IRS  sensor  errors  (i.e.  the  slowly  varying  errors  in  the  gyros  and  accelerometers)  show  a 
+/- 1 -sigma  bound  tracking  along  with  the  actual  filter  state  estimate  (rather  than  a  +/- 1 -sigma  bound  about 
zero)  to  give  an  idea  of  the  uncertainty  of  that  particular  estimate.  This  is  the  normal  graphical  display 
option  when  the  true  error  state  is  not  available  for  reference  purposes  (as  would  be  the  case,  normally,  for 
real  data). 


The  final  plots  to  be  shown,  based  on  this  Kalman  filter  run  (Figs.  24  and  25),  relate  to  the  so- 
called  residual  sequence  (or  innovations  sequence  -  see  Ref.  8)  for  each  of  the  5  measurement  time 
series  that  have  already  been  shown  in  Figs.  15  and  16.  For  a  given  scalar  measurement  time  series  Zj(k)  (i 
=  1,..,5),  the  innovations  time  series  process  associated  with  it  can  be  expressed  as 
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M<k)  =  Zj(k)  -  Hj(k)  5x'(k)  (20) 

If  the  filter  is  optimally  tuned,  then  Kalman  filter  theory  states  that  this  innovations  process  will  be  ZMWG 
and  will  have  the  theoretical  statistical  variance  computed  by  the  filter  as 

pta(k)}  =  Hi<k)P'(k)HjT(k)  +  r  j  (k)  (21) 


Figures  24  and  25  indeed  show  that  the  innovations  time  series  sequences  associated  with  each 
of  the  five  scalar  measurements  being  processed  by  the  filter  (i.e.  both  the  position  and  velocity-based 
measurements)  are  very  close  to  ZMWG  random  processes  and,  furthermore,  the  +/- 1 -sigma  bounds 
being  computed  (shown  in  dashed  lines  superimposed  on  the  same  plots  as  the  innovations)  appear  to 
closely  approximate  the  true  standard  deviation  in  each  case.  This  is  clear  evidence  that  the  Kalman  filter 
software  package  is  performing  properly.  For  the  case  of  real  navigation  data,  monitoring  the  innovations 
sequences  during  the  operation  of  the  Kalman  filter  can  be  an  effective  way  of  detecting  sudden  changes 
in  the  noise  characteristics  of  the  measurement  processes  that  might  be  due  to  a  sudden  fault  in  one  of 
the  measurement  sensors  (see  Refs.  16, 17). 

4.5  Conclusions  from  the  Simulation  Studies 

Prior  to  working  with  real  navigation  data  from  the  Twin  Otter,  a  great  many  Monte  Carlo  simulation 
runs  of  the  Kalman  filter  were  conducted  under  various  simulated  conditions,  primarily  to  check  out  the 
expected  robustness  and  accuracy  of  the  filter.  The  following  points  summarize  the  results  of  the  analysis 
of  the  simulation  runs: 

a)  For  this  particular  Kalman  filter  configuration,  simulation  results  indicate  that  there  appears  to  be  a 
small  improvement  to  be  gained  by  using  measurement  prefiltering  on  both  the  Loran-C  and 
Doppler-based  measurements,  as  theory  suggests.  However,  this  improvement  is  not 
considered  sufficient  to  warrant  using  prefiltering  in  a  real-time  application. 

b)  There  is  virtually  no  difference  in  filter  performance  when  comparing  results  with  and  without  the 
Doppler-based  measurements  being  included.  This  suggests  that,  in  the  presence  of  the 
reasonably  accurate  Loran-C  position-based  measurements,  the  rather  inaccurate  Doppler 
velocity  measurements  contribute  very  little  to  the  estimation  of  the  various  IRS  error  states. 

c)  One  robustness  check  involved  using  only  an  8  Hz  version  of  the  IRS  digital  data,  rather  than  the 
usual  16  Hz  version.  A  comparison  of  Kalman  filter  runs  for  these  two  different  IRS  data  rates 
confirmed  that  there  was  virtually  no  difference  in  filter  performance. 

d)  Another  filter  robustness  check  consisted  of  comparisons  of  filter  performance  for  update 
intervals  longer  than  the  usual  10  seconds.  A  Kalman  filter  update  interval  of  20  seconds 
provided  results  almost  identical  to  those  for  the  standard  10  second  interval.  For  a  30  second 
update  interval,  a  significant  degradation  in  filter  performance  was  apparent. 

e)  The  concept  of  using  so-called  Visual  On-Tops  (VOTs)  as  infrequent  measurements,  whenever 
they  occur,  was  emulated  in  the  filter  simulation  software  in  the  manner  described  in  Ref.  3. 
Simulation  results  confirm  that  VOT  processing  is  being  done  properly  within  the  filter  software, 
but  some  care  must  be  taken  as  to  the  details  of  how  VOT  measurements  are  processed  in  order 
to  avoid  filter  transients. 
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The  foregoing  robustness  checks  have  confirmed  the  inherent  soundness  of  the  Kalman  filter 
design  and  its  ability  to  degrade  gracefully  in  the  presence  of  computer  processing  shortcuts  that  might 
be  required  for  a  real-time  version  of  the  filter.  By  using  the  simulation  software  in  this  manner,  various 
tradeoffs  can  be  made  between  the  ultimate  in  filter  accuracy  (with  its  associated  computational 
complexity)  and  acceptable  accuracy  (with  significant  simplifications  in  the  filtering  algorithms). 

5.0  KALMAN  FILTER  RESULTS  USING  REAL  FLIGHT  DATA 

5.1  Specialized  Navigation  Flights 

Two  sets  of  specialized  navigation  flight  test  data  (each  approximately  1 .5  hours  in  duration)  were 
collected  onboard  the  Twin  Otter  for  use  in  evaluating  and  fine  tuning  the  Kalman  filter.  For  these  two 
flights,  a  route  was  chosen  such  that  the  aircraft  overflew  selected  landmarks  that  were  to  be  used  as 
reference  positions,  or  Visual  On-Tops  (VOT's).  This  provided  an  assessment  of  IRS  position  error  drift, 
Loran-C  navigational  accuracy,  and  Kalman  filter  performance  (in  regards  to  horizontal  positioning 
accuracy,  at  least)  for  the  entire  flight.  As  part  of  each  of  these  flights,  a  wind  box  and/or  straight  line  runs 
on  reciprocal  headings  along  a  railroad  line  were  included  in  order  to  evaluate  wind  computation  accuracy 
when  filter-corrected  IRS  velocity  components  were  being  used  in  the  wind  equations  (see  Ref.  18  for 
details  of  this  aspect  of  the  development  work).  Figure  26  shows  the  horizontal  track  plot  for  the  second 
of  these  flights  (Flight  #2)  -  the  numbers  from  1  to  15  correspond  to  the  locations  of  the  VOTs  along  the 
route  (some  VOTs  were  overflown  more  than  once  during  a  flight).  The  first  flight  (Flight  #1)  followed  a 
trajectory  identical  to  the  one  shown  in  Fig.  26,  except  that  the  wind  box  was  omitted.  Table  5.1  lists  the 
set  of  ten  individual  landmarks  that  were  used  as  VOTs,  together  with  their  geographical  locations  which 
are  felt  to  be  accurate  to  better  than  0.1  nm  RMS  in  each  of  their  latitude  and  longitude  components. 

TABLE  5.1 

GEOGRAPHICAL  LOCATIONS  OF  THE  VISUAL  ON-TOPS 


LANDMARK  DESCRIPTION 

LATITUDE 

LONGITUDE 

1) 

RUNWAY  INTERS.  -  OTTAWA  AIRPORT 

45.31667° 

75.66500° 

2) 

OTTAWA  VOR 

45.44167° 

75.89667° 

3) 

FITZROY  DAM 

45.47167° 

76.23333° 

4) 

CHENAUX  DAM 

45.58333° 

76.67500° 

5) 

MOUNT  ST.  PATRICK  TOWER 

45.31000° 

76.89167° 

6) 

MOUNTAIN  CHUTE  DAM 

45.19667° 

76.90833° 

7) 

BARRET  CHUTE  DAM 

45.24333° 

76.75833° 

8) 

BURNSTOWN  BRIDGE 

45.38500° 

75.57833° 

9) 

WEST  END  -  RAILROAD  RUN 

45.07667° 

75.90333° 

10) 

EAST  END  -  RAILROAD  RUN 

45.16000° 

75.83833° 
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5.2  Fundamental  Filter  Performance  Using  Real  Flight  Data 
5.2.1  Data  from  Flight  #1  (May  *88) 

Navigation  data  from  Flight  #1  were  the  first  data  from  the  Twin  Otter  to  be  analyzed,  in  detail,  from 
the  point  of  view  of  suitability  for  use  in  the  Kalman  filter  algorithms  that  had  been  developed.  Initial 
analysis  and  use  of  these  data  in  the  Kalman  filter  program  revealed  a  few  minor  problems  that  had  to  be 
overcome.  The  plots  of  the  two  position  measurement  components  (i.e.  and  z2,  based  on  the  Loran-C 
position  data  being  differenced  with  the  IRS  position  data),  shown  in  Fig.  27,  indicate  two  of  these 
problems.  A  very  obvious  'spiking'  phenomenon  occurs  in  the  time  series  plot  of  the  longitude 
measurement  component.  This  effect  was  traced  to  a  fundamental  flaw  in  the  onboard  software  for 
formatting  the  data  to  be  recorded  that  caused  an  intermediate-level  bit  in  the  IRS  longitude  digital  word  to 
stick  at  certain  times  -  this  software  bug  was  fixed  for  subsequent  flights.  In  the  case  of  this  corrupted  set 
of  data,  a  simple  on-line  algorithm  was  used  in  order  to  detect  and  smooth  out  the  spikes  prior  to  the  data 
being  used  in  the  Kalman  filter.  Notice,  also,  that  in  the  time  series  plots  for  each  of  the  measurement 
components  a  rather  sudden  shift  of  about  0.3  nm  occurs  several  times  throughout  the  flight.  A  thorough 
analysis  of  this  phenomenon  revealed  that  the  measurement  shifts  always  occurred  during  times  when 
the  aircraft  was  turning.  At  first  it  was  thought  that  the  Loran-C  receiver  might  be  sensitive  to  aircraft 
manoeuvres  and  build  up  an  error  during  turns.  However,  experience  using  this  particular  ARNAV  Loran- 
C  receiver  onboard  both  the  Twin  Otter  and  the  Convair  580  has  revealed  the  existence  of  a  consistent  4 
second  time  lag  in  the  output  from  the  receiver  that  becomes  particularly  noticeable  during  manoeuvres. 
This  was  corroborated  by  time-shifting  the  Loran-C  data,  relative  to  the  IRS  data,  by  exactly  4  seconds  -  the 
spurious  shifts  in  the  position  measurement  data  disappeared.  Figure  28  shows  the  same  measurement 
data  after  applying  the  corrections  that  have  been  described  above  -  each  of  the  time  series  traces  is  now 
relatively  smooth,  with  only  the  near-constant  level  of  high  frequency  measurement  noise  (estimated  to 
be  about  0.07  nm  RMS  for  each  component)  now  remaining.  Figure  29  shows  the  corresponding  time 
series  traces  for  the  three  velocity  measurements  (i.e.  z3  to  Z5,  based  on  the  Doppler  velocity  sensing  data 
being  differenced  with  the  IRS  velocity  data).  As  can  be  seen  in  Fig.  29,  Doppler-based  measurements 
are  quite  noisy,  due  to  the  rather  high  levels  of  fluctuation  no.oe  inherent  to  this  particular  Doppler  velocity 
sensor.  Moreover,  the  levels  of  the  fluctuation  noise  are  somewhat  different  for  each  of  the  U,  V  and  W 
velocity  components  -  i.e.  about  3.0  ft/s  RMS  for  U,  10.0  ft/s  RMS  for  V  and  5.0  ft/s  RMS  for  W. 

The  ten  VOT's  listed  in  Table  5.1  were  used  to  determine  the  long-term  positioning  accuracy  of 
both  the  LTN-90-100  IRS  and  Loran-C  receiver  during  the  course  of  the  flight.  Figure  30  is  a  set  of  plots 
depicting  the  buildup  in  IRS  positional  error  (latitude  and  longitude  components)  during  the  course  of 
Flight  #1 ,  as  determined  by  position  comparisons  at  the  visual  landmark  locations.  Positional  error  drifts  of 
about  1  nm/hour  are  typical  for  this  quality  of  IRS,  together  with  a  Schuler-induced  oscillating  error  (with  an 
84  minute  period)  of  perhaps  1  nm  peak-to-peak.  There  is  evidence  of  a  very  significant  Schuler  error 
component  building  up  in  the  latitude  error  trace  of  Fig.  30,  for  example.  In  contrast,  Fig.  31  shows  a 
similar  series  of  plots  for  the  Loran-C  position  errors  -  note  the  small,  bounded  nature  of  these  errors,  with 
both  components  consistently  below  0.2  nm  in  magnitude.  These  error  levels  are  consistent  with  Loran-C 
system  specifications  for  good  coverage  areas,  such  as  the  area  where  these  data  were  collected.  The 
small,  bounded  nature  of  Loran-C  errors  makes  airborne  Loran-C  an  ideal  redundant  navaid  for  identifying 
the  dominant  IRS  errors  using  a  Kalman  filter  integrated  navigation  scheme. 

Figure  32  is  a  listing  of  the  SDKFNAV  CTRL  file  used  for  controlling  the  running  of  the 
SDKFNAV  software  package  on  this  particular  set  of  flight  data.  The  statistical  specifications  in  this  file 
have  been  fine  tuned  for  the  known  noise  levels  of  each  of  the  measurement  types  and  for  a  slightly  out- 
of-specification  situation  regarding  the  x  accelerometer  of  the  IRS  (identified  by  the  manufacturer  during  a 
subsequent  re-calibration  exercise).  Also,  the  VOT  lat/long  position  values  have  been  time  shifted  so  that 
they  are  precisely  synchronized  with  the  10  second  update  intervals  of  the  Kalman  filter.  Figures  33  to  39 
show  the  time  series  plots  of  the  complete  set  of  error  states  estimated  by  the  Kalman  filter  based  on  the 
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control  conditions  specified  in  SDKFNAV  CTRL  (i.e.  Fig.  32).  For  the  filter  run  depicted  in  these  plots, 
the  specified  measurement  noise  levels  for  each  of  the  Loran-based  and  Doppler-based  measurement 
types  were  fine-tuned  to  produce  the  best  overall  results.  These  plots  also  include  dashed  lines  that 
indicate  the  ± to  bounds  computed  by  the  Kalman  filter  to  reflect  the  accuracy  level  of  each  error  state 
estimate.  Note  that  the  error  states  associated  with  the  vertical  channel  (i.e.  IRS  altitude,  IRS  vertical 
velocity,  IRS  vertical  acceleration  and  barometric  altimeter  bias)  remain  very  close  to  zero  because  these 
errors  cannot  really  be  observed  by  the  Kalman  filter  in  the  absence  of  an  accurate,  independent  source 
of  altitude  information  to  use  as  a  redundant  measurement  type.  As  well,  the  sea  bias  components  are 
zeroed  out  because  no  overwater  flying  occurred  during  the  course  of  Flight  #1 .  Figure  40  is  a  plot  of  the 
Kalman  filter's  IRS  position  calibration  accuracy  at  each  of  the  VOT’s  occurring  during  the  flight,  together 
with  ±lo  bounds  to  describe  the  accuracy  predicted  by  the  filter.  The  uncertainty  bounds  settle  out  at 
approximately  ±  0.2  nm,  and  the  true  error  of  the  estimates  always  falls  within  the  bounds  -  a  good 
indication  of  a  robust  Kalman  filter  design.  Time  series  plots  of  the  innovations  sequences  corresponding 
to  each  of  the  5  scalar  measurements  (shown  in  Figs.  40a  and  40b)  appear  to  be  close  to  ZMWG  random 
sequences,  which  is  another  confirmation  of  proper  Kalman  filter  performance.  Although  there  is  no 
independent,  accurate  measure  of  aircraft  inertial  velocity  to  use  in  order  to  verify  the  filter’s  estimation 
accuracy  for  IRS  velocity  errors,  it  should  be  noted  that  the  filter  estimates  of  IRS  velocity  errors  are 
reasonably  consistent  when  compared  to  the  slope  of  the  position  error  traces  of  Fig.  23  (these  traces  are 
known  to  be  accurate  to  within  ±0.2  nm).  All  of  the  other  error  states  show  consistent  behaviour  and  attain 
values  in  line  with  those  predicted  during  simulation  studies. 

5.2.2  Data  from  Flight  #2  (May  '89) 

Figure  41  shows  the  time  series  plots  of  the  2  Kalman  filter  position  measurement  components 
(i.e.  zi  and  z2)  that  are  based  on  the  navigation  data  from  Flight  #2.  Note  that,  as  in  the  case  with  Flight  #1 , 
the  Loran-C  data  has  been  time-shifted  by  4  seconds  in  order  to  synchronize  it  with  the  IRS  data.  The 
near-constant  high  frequency  measurement  noise  levels  are  about  the  same  as  for  the  corresponding 
data  from  Flight  #1  (see  Fig.  28),  namely  0.07  nm  RMS  for  each  of  the  components.  Figure  42  shows  the 
time  series  plots  for  the  3  Kalman  filter  velocity  measurement  components  (i.e.  Z3to  Z5)  -  the  U,  V,  and  W 
Doppler-based  velocity  measurements  look  much  the  same  as  the  corresponding  plots  for  the  Flight  #1 
data  (i.e.  Fig.  29)  and  have  the  same  levels  of  high  frequency  measurement  noise  as  on  the  earlier  Flight 
#1 .  With  regard  to  the  Doppler  system,  a  very  minor  problem  occurred  with  it  during  the  flight  -  it  went  into 
memory  for  about  20  seconds  early  on,  due  to  flying  over  smooth  water  for  a  brief  period  of  time.  The 
simple  solution  to  this  problem  was  to  incorporate  logic  in  the  Kalman  filter  program  to  skip  the  processing 
of  Doppler-based  measurement  data  around  the  point  in  time  when  the  Doppler  data  were  suspect. 
Figure  43  shows  the  IRS  horizontal  position  errors  at  the  VOT  reference  positions  during  the  flight.  IRS 
position  error  drifting  was  not  quite  as  bad  as  what  had  occurred  during  Flight  #1 ,  due  to  the  fact  that  an 
out-of-spec  accelerometer  had  now  been  repaired.  For  the  same  set  of  VOT  reference  landmarks,  Fig.  44 
shows  the  Loran-C  position  errors  during  the  course  of  the  flight ,  with  results  that  are  much  the  same  as 
those  for  Flight  #1  -  i.e.  the  errors  in  the  Loran-C  position  components  are  seen  to  be  bounded  in 
magnitude  by  0.2  nm. 

Figures  45  to  51  show  the  complete  set  of  error  state  time  histories  and  associated  ±1a  bounds 
for  a  Kalman  filter  run  using  the  Flight  #2  data  and  an  SDNAV  CTRL  file  almost  identical  to  the  one 
shown  in  Fig.  32  -  the  only  difference  being  that  the  statistics  for  the  x  accelerometer  have  been  changed 
to  reflect  the  fact  that  this  accelerometer  has  been  re  calibrated  and  is  back  within  specification.  The 
overall  results  are  quite  similar  to  those  for  the  corresponding  filter  run  on  the  data  from  Flight  #1  (i.e.  Figs. 
33  to  39)  apart  from  the  obvious  fact  that  IRS  horizontal  position  and  velocity  errors  are  now  a  lot  smaller. 
Finally,  Fig.  52  shows  the  filter's  IRS  horizontal  position  calibration  accuracy  at  each  of  the  VOT  landmark 
points  during  the  flight,  together  with  the  associated  ±lo  uncertainty  bounds  that  confirm  that  actual  error 
levels  always  fall  within  the  bounds,  as  desired. 
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5.3  A  Comparison  of  Three  Different  Measurement  Scenarios 

There  are  three  different  measurement  scenarios  that  can  be  analyzed  easily,  based  on  different 
combinations  of  the  navaid  data  being  used  as  Kalman  filter  measurements.  The  three  configurations  to 
be  considered  are  the  following:  i)  only  the  Loran-C  position  measurements  available,  ii)  both  the  Loran-C 
and  Doppler-based  measurements  available  and  iii)  only  the  Doppler-based  velocity  measurements 
available.  Kalman  filter  IRS  position  error  estimation  results  for  each  of  the  three  filter  configurations,  and 
both  sets  of  navigation  data,  are  shown  in  Figs.  53  to  56.  A  comparison  of  position  error  estimation  for 
each  of  the  three  filter  configurations  reveals  that  case  i)  and  case  ii)  results  are  almost  identical,  with 
position  errors  consistently  below  0.2  nm  in  magnitude  (as  indicated  by  the  uncertainty  bounds),  while 
case  iii)  results  (i.e.  Doppler-only)  degrade  to  about  0.6  nm  RMS.  Kalman  filter  IRS  velocity  error  estimation 
results  for  each  of  the  three  filter  configurations,  and  both  sets  of  flight  data,  are  shown  in  Figs.  57  to  60. 
The  case  i)  results  (Loran-C  only)  are  the  most  accurate,  due  to  the  high  quality  nature  of  the  Loran-C 
measurement  data.  Next  in  accuracy  would  be  the  case  ii)  results  (Loran-C/Doppler  measurement 
combination)  -  there  is  a  bit  of  accuracy  degradation  due  to  the  fact  that  the  inferior  quality  of  the  Doppler 
velocity  data  makes  accurate  Doppler  error  modelling  somewhat  difficult  and  slightly  corrupts  the 
estimation  results,  but  the  error  traces  are  fundamentally  the  same  as  in  case  i).  Also,  in  both  of  these 
cases,  the  RMS  error  bounds  settle  out  at  about  1 .0  ft/s  (0.3  m/s)  for  each  component.  Finally,  the  case  iii) 
(Doppler  only)  results  are  the  least  accurate,  as  judged  by  the  wider  RMS  bounds  (±  2.3  ft/s  or  ±  0.7  m/s 
on  average)  and  sensitivity  to  aircraft  manoeuvres.  However,  even  in  the  Doppler-only  measurement 
case,  the  dominant  trends  in  the  IRS  velocity  error  traces  of  case  i)  are  retained.  The  foregoing 
comparisons  serve  to  point  out  that,  even  in  the  absence  of  the  Loran-C  position  data,  there  is  still  a 
significant  improvement  in  inertial  position  and  velocity  accuracy  to  be  gained  from  the  implementation  of 
this  Kalman  filter  design  using  Doppler-only  measurements. 

5.4  Conclusions  from  the  Analysis  of  Real  Flight  Data 

1 )  Loran-C  coverage  in  the  vicinity  of  Ottawa  (i.e.  using  the  NE  USA  chain)  is  excellent,  with  actual 
Loran-C  position  errors  well  within  the  stated  Loran-C  performance  specification  of  0.25  nm  RMS 
for  each  component.  As  a  result,  with  Loran-C  measurements  the  Kalman  filter  can  consistently 
provide  IRS-based  positioning  accuracy  to  a  level  of  0.2  nm  RMS  (i.e.  the  same  accuracy  level  as 
Loran-C)  and  IRS-based  velocity  accuracy  to  a  level  of  1 .0  ft/s  RMS  in  the  horizontal  components. 
IRS  attitude  errors  can  also  be  decreased  somewhat,  from  0.05  deg  RMS  down  to  0.03  deg  RMS 
in  pitch/roll  and  from  0.2  deg  RMS  down  to  0.08  deg  RMS  in  heading. 

2)  With  the  vertical  channel  of  the  IRS  being  controlled  by  an  internal,  third-order  baro-damping  loop, 
vertical  position  accuracy  is  determined  by  the  fundamental  barometric  altimeter  accuracy,  and  IRS 
vertical  velocity  accuracy  is  better  than  0.5  ft/s  RMS  without  any  Kalman  filtering  being  applied.  In 
the  absence  of  a  redundant  source  of  altitude  information,  substantially  more  accurate  than 
barometric  altitude,  there  is  no  further  improvement  to  IRS  vertical  channel  accuracy  that  can  be 
provided  by  the  present  version  of  the  Kalman  filter.  However,  the  present  Kalman  filter  design  is 
readily  adaptable  if  and  when  redundant  vertical  position  and/or  velocity  measurements  become 
available  (as  would  be  the  case  if  GPS  were  included  in  the  suite  of  navigation  sensors). 

3)  Initial  attempts  at  applying  Kalman-corrected  IRS  inertial  velocity  data  to  the  wind  computation  task 
have  proved  encouraging  (see  Ref.  18  for  more  details).  Wind  calculations  using  corrected  IRS 
velocities  appear  to  give  more  consistent  results,  and  wind  computation  accuracy  of  1  ft/s  RMS 
per  channel  should  be  attainable. 
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4)  It  has  been  established  that  there  is  no  lurther  improvement  in  IRS  error  estimation  accuracy  to  be 
gained  by  processing  the  Doppler  measurements  along  with  the  Loran-C  measurements.  As  a 
matter  of  fact,  with  this  particular  Doppler  velocity  sensor,  one  must  be  careful  to  establish  the 
measurement  noise  'weighting'  property  so  as  to  avoid  having  the  Doppler-based  measurement 
data  actually  degrade  filter  performance  when  it  is  used  in  conjunction  with  the  Loran-C  data.  A 
Kalman  filter  analysis  of  the  two  sets  of  specialized  navigation  flight  test  data  has  revealed  that 
Doppler  U  channel  scale  factor  errors  can  be  1 .5%  -->  2%,  even  after  careful  calibration  of  the 
system,  and  are  quite  variable  because  of  the  terrain  sensitivity  of  this  rather  old  Doppler  set.  The 
V  channel  boresight  error  is  not  particularly  constant  either,  with  variations  that  that  can  exceed  1 
deg  at  times.  With  the  200  ft/s  groundspeed  typically  attained  by  the  Twin  Otter,  uncalibrated 
Doppler  velocity  errors  in  each  of  the  U  and  V  channels  can  be  3.0  ft/s  or  more.  However,  it  is 
worthwhile  to  retain  Doppler  velocity  measurements  in  the  Kalman  filter  configuration  for  situations 
in  which  Loran-C  coverage  is  intermittent  or  simply  not  available.  For  the  Doppler-only 
measurement  case,  IRS  horizontal  position  accuracy  is  approximately  0.6  nm  RMS  and  horizontal 
velocity  accuracy  is  about  2.3  ft/s  RMS.  Due  to  the  error  levels  inherent  to  the  Doppler  data,  IRS 
attitude  errors  cannot  be  reduced  by  any  significant  amount  when  only  Doppler  data  are  available. 

5)  The  IRS  gyro  and  accelerometer  bias-like  error  states  require  a  much  longer  time  period  than  1 .5 
hours  to  be  estimated  accurately  by  the  Kalman  fitter  -  i.e.  more  like  3  to  4  hours  (the  equivalent  of 
a  few  periods  of  the  84  minute  Schuler  oscillation).  However,  even  in  the  short  term,  IRS  sensor 
error  state  estimation  does  appear  to  be  yielding  consistent  intermediate  values,  judging  by  the 
nature  of  the  time  histories  for  these  errors  over  the  first  1 .5  hours. 

6)  Analysis  of  various  Kalman  fitter  runs,  using  a  measurement  prefiltering  option  for  both  the  Loran- 
C  and  Doppler  data,  has  revealed  that  estimation  accuracy  does  not  improve  significantly.  It  is 
concluded  that  measurement  prefiltering  is  not  worth  the  extra  computational  burden  for  on-line 
filtering  applications,  but  should  be  retained  as  an  option  for  off-line  Kalman  filter  analysis.  Various 
results  using  simulated  data  support  this  same  conclusion. 

7)  The  Kalman  filter  performed  almost  identically  when  using  an  8  Hz  version  of  the  IRS  data,  rather 
than  the  usual  16  Hz  version.  This  is  encouraging  from  the  standpoint  of  designing  a 
computationally  efficient  version  of  the  filter  for  real-time  operations.  Further  experimentation 
revealed  that  it  is  quite  feasible  to  'slow  down’  the  filter  to  a  0.05  Hz  update  rate  (i.e.  20  second 
update  interval)  without  a  significant  degradation  in  performance. 

8)  Linder  certain  circumstances,  using  VOT's  as  extra  measurements  (whenever  they  occur)  can 
cause  an  undesirable  transient  in  several  of  the  Kalman  filter  error  states.  The  process  of  utilizing 
VOT's  as  measurements  should  be  checked  out  very  carefully,  for  any  given  Kalman  filter 
configuration,  to  ensure  that  basic  filtering  integrity  is  maintained. 


6.0  CONCLUSIONS  AND  FUTURE  WORK 

The  Kalman  filter  integrated  navigation  design,  described  in  detail  in  this  report,  has  been 
demonstrated  to  meet  the  performance  specifications,  when  using  real  flight  data,  that  had  been 
determined  from  simulation  studies.  Most  importantly,  filter-corrected  IRS  position  and  velocity 
components  have  been  verified  to  be  consistently  within  their  predicted  accuracy  limits.  In  particular,  for 
the  Twin  Otter  wind  computation  task,  a  guaranteed  IRS  velocity  error  limit  of  1  ft/s  will  be  a  very  significant 
improvement  in  performance  compared  to  the  present  error  limits  of  3  to  4  ft/s  for  Doppler  velocity  and  3  to 
6  ft/s  for  raw  IRS  velocity.  Applying  this  Kalman  filter  approach  to  airborne  wind  measurement  will  ensure 
that  the  dominant  source  of  error  will  no  longer  be  the  sensing  of  inertial  velocity  but,  rather,  the 
measurement  of  the  air  data  information.  Any  inconsistencies  and  errors  in  wind  data  that  show  up  during 
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analysis  can  then  be  attributed  solely  to  the  air  data  system  and  appropriate  action  taken.  The  Kalman  filter 
software  described  in  this  report  has  been  successfully  ported  to  the  VAX/VMS  environment.  Recent 
flight  data  (i.e.  from  the  summer  of  1990),  acquired  using  the  latest  version  of  the  Twin  Otter's  data 
acquisition  software,  has  been  validated  to  be  completely  compatible  with  the  VMS  Kalman  filter  software 
package.  The  capability  of  providing  routine  Kalman  filter-based  IRS  corrections  as  a  standard  feature  of 
the  Twin  Otter  playback  software  is  ready  to  be  tested  out  on  the  latest  data  from  Twin  Otter  field 
experiments. 

Future  work  will  involve  confirming  the  accuracy  and  reliability  of  the  IRS  correction  software,  using 
various  sets  of  navigation  data  from  Twin  Otter  field  experiments,  by  clearly  demonstrating  consistent 
improvement  in  the  wind  computations  that  are  based  on  Kalman  filter-corrected  IRS  velocity  data. 
Software  will  be  developed  for  a  Kalman  filter/smoother  algorithm  that  is,  potentially,  more  accurate  than 
the  present  Kalman  filter  algorithm  for  post-processing  Twin  Otter  data.  A  more  powerful  airborne 
microprocessor  (a  DEC  MicroVax-ll)  will  soon  be  installed  onboard  the  Twin  Otter.  Once  that  occurs,  the 
Kalman  filter  software  will  be  installed  in  it  and  employed  for  real-time  wind  compuations.  As  it  so  happens, 
the  present  version  of  the  ground-based  Kalman  filter  software  can  be  easily  employed  for  the  real-time 
task.  At  some  point  in  the  near  future,  the  feasibility  of  installing  a  GPS  receiver  onboard  the  Twin  Otter  will 
be  assessed  with  a  view  to  gaining  some  experience  with  an  IRS/GPS  Kalman  filter  configuration.  This 
sensor  combination  has  the  potential  for  providing  inertial  data,  at  a  16  Hz  rate,  with  a  real-time  absolute 
accuracy  of  30  feet  in  each  position  component,  0.3  ft/s  in  each  velocity  component,  0.005  deg  in  each  of 
pitch/roll  attitude  and  0.02  deg  in  heading  (see  Refs.  19,  20).  For  certain  contingency  situations  that 
might  occur  while  running  a  Kalman  filter,  a  careful  study  will  be  conducted  into  the  best  way  to  handle 
unwanted  transient  errors  in  a  practical  manner.  Two  common  contingency  situations  in  which  this 
transient  phenomenon  can  manifest  itself  are  i)  the  occurrence  of  a  VOT  landmark  along  the  flight  path,  to 
be  processed  as  an  auxiliary  position  measurement  and  ii)  the  case  where  one  of  the  redundant  sources 
of  filter  measurement  data  drops  oft  line  at  some  point  during  the  filter  run.  One  solution  to  this  filter  re¬ 
configuration  problem  that  is  being  looked  at  very  seriously  is  to  run  several  Kalman  filters  in  parallel,  each 
one  depicting  a  possible  Kalman  filter  configuration  based  on  the  various  sources  of  navigation  data  that 
might  be  available.  The  outputs  from  this  bank  of  Kalman  filters  would  be  blended  together,  using 
weightings  that  would  depend  on  real-time  information  about  the  integrity  of  the  various  navigation 
sensors  that  comprise  the  complete  navigation  suite. 
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APPENDIX  A 


Fundamental  Inertial  Navigation  Definitions  and  Equations 


The  basic  inertial  system  parameters  and  equations  are  defined  in  this  Appendix  as  background 
material  for  the  details  of  INS/IRS  error  modelling  and  associated  Kalman  filter  design  that  follow  in 
subsequent  Appendices. 


a)  Model  of  Earth’s  Geoid  - 

The  earth's  reference  ellipsoid  (Ref.  14)  can  be  defined  in  terms  of  the  following  three 
parameters: 


R0  (equatorial  radius):  3441 .726  nm,  6378.165  km  or  2.0925732  x  107  ft 
e  (earth's  flattening  or  ellipticity):  1/298.25  =  0.0033529 
e2  (earth's  eccentricity):  e2  =  2e-e2  =  0.00669460 


Define  L  to  be  geographic  latitude  and  X  to  be  geographic  longitude,  with  north  latitude  and  west 
longitude  being  the  positive  senses  (note  that,  usually,  east  longitude  is  considered  positive).  Two 
important  radii  of  curvature  are  of  interest,  namely: 


q.  (sometimes  denoted  pm):  the  meridian  radius  of  curvature 
r\  (sometimes  denoted  pp):  the  prime  radius  of  curvature 


Approximate  expressions  for  these  radii,  suitable  for  INS  work,  are  given  as: 


H.  =  Ro  [  1  -  e2  +  1 .5  e2  sin2  L  ] 

f\  =  RoM  +  (e2/2)  sin2  L )  (A1) 


b)  Model  of  Earth’s  Gravity  - 

Let  h  represent  altitude  (in  feet)  above  the  reference  ellipsoid.  An  adequate  mathematical  model 
for  so-called  plumb  bob  gravity  go  (i.e.  including  earth's  rotation  effects)  is  given  by 

go  -  Gi(  1  - 2h/R0  +  2e sin2  L )  +  G2(1-3sin2L)  -  G3(1-esin2L  +  h/R0)cos2L  (A2) 

with  Gi  =  32. 14647 ft/s2,  G2  =  0.0521882  ft/S2  and  G3  =  0.1106653  ft/S2. 
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c)  Velocity  and  Position  Equations  In  Geographic  Coordinates  - 

We  want  to  express  aircraft  velocity  and  position  components,  relative  to  the  rotating  earth,  using 
the  local-level  geographic  coordinate  frame  (i.e.  N,  E  and  z)  with  z  positive  upwards.  First  of  all,  several 
vectors  expressed  in  the  geographic  coordinate  frame  must  be  defined: 

v(vN,vE,vz)  -  aircraft  velocity  relative  to  the  earth 

P(Pn.  Pe.  Pz)  *  angular  velocity  of  geographic  frame  relative  to  the  earth 

f  (fN,  f e  .  fz)  -  specific  force  vector  (i.e.  net  acceleration  other  than  that  due  to  gravity) 

9  (9n.  9e.  9z)  -  plumb  bob  gravity  vector 

Q  (Qn  ,  QE  .  Qz)  -  angular  velocity  of  the  earth  relative  to  inertial  space 
co((on,  coe  ,  ®z)  -  angular  velocity  of  geographic  frame  relative  to  inertial  space 

Let  |Q|  denote  the  magnitude  of  the  earth  rate  vector  Q  (the  magnitude  is  7.2921  x  10  5  rad/s).  The 
components  of  Q  are  then 

On  =  |Q(cosL;  Qe  =  0;  =  |Q|sinL  (A3) 

Let  p  be  the  angular  velocity  of  the  geographic  frame  relative  to  the  earth  which  can  be  expressed  in 
component  form  as 

PN  =  vE/(rx  +  h);  pe  =  -vN/(rL  +  h);  pz  =  vE  tan  L7(rx  +  h)  (A4) 


Then,  the  angular  velocity  of  the  geographic  frame  relative  to  inertial  space  to  can  be  expressed  in 
component  form  as 


c^j  =  Pn  +  On:  <»e  =  Pe:  ®z  =  Pz  + 


(A5) 


It  is  also  useful  to  note  the  following  expansions: 

(flN  +  On  =  pN  +  2  Qn  =  pN  +  2  |Q|  COS  L 

coz  +  Qz  =  Pz  +  2  Qz  =  pz  +  2  )A|  sin  L  (A6) 


With  the  foregoing  definitions,  the  differential  equations  describing  inertial  velocity  dynamics  are  given  by 
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VN  =  *N  +  9N  -  (a>z  +  Qz)VE  +  (0£VZ 
VE  =  +  9e  -  ((On+An )Vz  +  ((Oz+nzJVN 

vz  =  fz  +  9z  -  vN  +  (con  +  £2n)  ve 


(A7) 


and  it  is  normally  assumed  that  gN  =  gE  =  0,  gz  =  -gD-  The  inertial  latitude,  longitude  and  height  dynamics 
are  governed  by  the  following  differential  equations: 


L  =  vN/(rL  +  h) 

X  =  -V£/[(rx  +h)cosL] 

h  =  vz  (A8) 


d)  Baro-lnertlal  Vertical  Channel  - 

It  is  well  known  that  the  vertical  navigation  equations  (i.e.  vz  in  Eqn.  A7  and  h  in  Eqn.  A8)  are 
inherently  unstable  due  to  the  particular  form  of  gz  =  -go(h,L).  A  small  positive  error  in  h  results  in  a  small 
negative  error  in  go  and,  hence,  a  small  positive  error  in  gz.  This  leads  to  a  positive  error  in  vz,  which 
integrates  into  a  larger  and  larger  error  in  h  -  the  instability  is  exponential  with  a  time  constant  of  about  ten 
minutes.  This  instability  is  damped  out  in  the  LTN-90-100  IRS  using  a  baro  altimeter  input  hb  in  a  third- 
order  digital  filtering  loop  (see  Ref. 9).  The  modified  versions  of  the  h  and  vz  equations  to  define  a  third- 
order  loop  are 


h  =  vz  -  Ki  (h-l\) 

vz  =  fz  +  gz  -  «e  vn  +  (con  +  *2n)  ve  -  K2  (h  -  hb)  -  5a 

5a  =  K3(h-hb)  (A9) 

where  5a  is  a  vertical  acceleration  correction  term.  A  triple  pole  placement  at  s  =  -1/x  (x  is  the  time  constant 
of  the  loop)  occurs  in  the  third-order  baro-inertial  channel  if 

K,  =  3/x;  K2  =  3/x2  +  2g/R0 ;  K3  =  1/x2  (A10) 


where  g  is  the  nominal  value  of  earth's  gravity.  In  Litton's  LTN-90-100  IRS,  x  has  been  set  to  20  seconds, 
resulting  in  coefficient  values  of  Ki  =  0.15,  K2  =  0.75  x  10-2  and  K3  =  0.125  x  10  3. 
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e)  Direction  Cosine  Matrix  (PCM)  - 

An  important  coordinate  transformation  for  the  strapdown  INS  case  is  the  one  from  aircraft  body 
axes  (x,  y,  z)  to  local-level  geographic  coordinates  (N,  E,  z)  via  the  so-called  Direction  Cosine  Matrix  (DCM) 
transformation  CbG.  Aircraft  body  axes  are  defined  as  i)  x:  longitudinal  (forward  positive),  ii)  y:  lateral  (right 
positive)  and  iii)  z:  normal  (down  positive).  Local-level  geographic  coordinates  are  defined  as  i)  N:  True 
North  positive,  ii)  E:  True  East  positive  and  iii)  z:  local  vertical  (up  positive).  The  elements  of  CbG  can  be 
expressed  as  transcendental  functions  of  the  Euler  angles  (6, 4>  and  y)  as  follows: 

Ci  i  =  oos0  cosy 

Ci  2  =  sin 0  sin 4> cosy  -  cos<)>  siny 
Ci 3  =  sin 0  cos 4»  cos y  +  sin  4>  siny 
C2,i  =  cos0  siny 

C22  =  sin0  sin 4> siny  +  cos <t>  cosy 
C2,3  =  sin  0  cos  4>  sin  y  -  sin  $  cos  y 
C31  =  sin0 
C3  2  =  -  cos  0  sin  $ 

C33  «  -COS 9  OOS(|) 


(All) 
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APPENDIX  B 


Error  State  Equations  for  the  IRS/Doppler/Loran-C  Kalman  Filter 


Details  of  the  error  model  for  the  complete  24-state  Kalman  filter  design,  in  terms  of  continuous¬ 
time  differential  equations,  are  given  in  this  Appendix.  By  explicitly  writing  out  the  error  state  equations  in 
the  desired  order,  the  individual  elements  of  matrices  F  and  G  can  be  readily  identified. 


a)  Continuous-Time  Equations  - 


The  error  state  vector/matrix  differential  equation  has  the  following  general  form: 


8x 


F  8x  +  G  u 


(B1) 


where,  in  the  case  of  this  particular  Kalman  filter  design,  the  error  state  vector  5x  has  the  24  elements 

8x  =  [8L  8k  8h  5vn  8ve  8vzeN  eg  e2  8a  Bco*  Ba>y  Btoz  Bax  Bay  Baz  Bhb  Blat  Blnq 
SFu  By  Byy  SBfl  SBg  ]  T 

(B2) 

and  the  ZMWG  plant  noise  vector  u  has  the  24  elements 

u  =  [  u(Dx  uo)y  uroj.  uax  uay  uaz  ugu  uge  ugz  uhb  uBcox  uBcoy  uBcoz  uBax  uBay  uBaz 

uBhb  uBlj  uBlg  uSFj  uBv  uBw  uSBn  uSBe]t  (B3) 


b)  Definition  of  Error  States  - 

A  very  short  description  of  each  of  the  error  states  is  given  below: 

1.  8L  -  error  in  IRS  latitude  (rad) 

2.  8A.  -  error  in  IRS  longitude  (rad) 

3.  8h  -  error  in  IRS  altitude  (feet) 

4.  Svn  -  error  in  :RS  north  velocity  component  (ft/sec) 

5 .  SvE  -  error  in  IRS  east  velocity  component  (ft/sec) 

6.  8vz  -  error  in  IRS  local  vertical  velocity  component  (ft/sec) 

7.  en  -  IRS  attitude  error  about  N  (north)  axis  (rad) 

8.  ee  -  IRS  attitude  error  about  E  (east)  axis  (rad) 

9 .  ez  -  IRS  attitude  error  about  z  (local  vertical)  axis  (rad) 
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10.  5a  -  acceleration  correction  in  baro  damped  vertical  channel 

1 1 .  B<0x  *  x  axis  strapdown  gyro  Markov  error  (rad/sec) 

12.  Bcoy  -  y  axis  strapdown  gyro  Markov  error  (rad/sec) 

13.  Btoz  -  z  axis  strapdown  gyro  Markov  error  (rad/sec) 

14.  Bax  -  x  axis  strapdown  accelerometer  Markov  error  (ft/sec2) 

15.  Bay  -  y  axis  strapdown  accelerometer  Markov  error  (ft/sec2) 

16.  Baz  -  z  axis  strapdown  accelerometer  Markov  error  (ft/sec2) 

17.  Bhb  -  baro  altimeter  Markov  error  (feet) 

18.  Blat  -  Loran-C  latitude  Markov  error  (rad) 

19.  Blng  -  Loran-C  longitude  Markov  error  (rad) 

20.  SFu  -  Doppler  U  channel  scale  factor  Markov  error  ( ) 

21.  By  -  Doppler  V  channel  boresight  Markov  error  (rad) 

22.  Bw  -  Doppler  W  channel  boresight  Markov  error  (rad) 

23.  SBn  -  North  Markov  sea  bias  component  (ft/sec) 

24.  SBe  -  East  Markov  sea  bias  component  (ft/sec) 


c)  Definition  of  Plant  Noise  Components  - 

The  24  elements  of  the  plant  noise  vector  u  can  be  described  as  follows: 

1.  uo)x  -  x  axis  gyro  random  drift  (rad/sec) 

2.  utoy  -  y  axis  gyro  random  drift  (rad/sec) 

3.  ucoz  -  z  axis  gyro  random  drift  (rad/sec) 

4.  uax  -  x  axis  accelerometer  random  noise  (ft/sec2) 

5.  uay  -  y  axis  accelerometer  random  noise  (ft/sec2) 

6.  uaz  -  z  axis  accelerometer  random  noise  (ft/sec2) 

7.  ugN  -  North  component  of  random  gravity  anomaly  (ft/sec2) 

8.  ugE  -  East  component  of  random  gravity  anomaly  (ft/sec2) 

9.  ugz  -  vertical  component  of  random  gravity  anomaly  (ft/sec2) 

10.  uhb  -  baro  altimeter  random  noise  (feet) 

1 1 .  uBrox  -  x  axis  gyro  Markov  driving  noise  (rad/s3'2) 

12.  uBtOy  -  y  axis  gyro  Markov  driving  noise  (rad/s3'2) 

1 3 .  uBtOz  -  z  axis  gyro  Markov  driving  noise  (rad/s3/2) 

14.  uBax  -  x  accelerometer  Markov  driving  noise  (ft/s5'2) 

15.  uBay  -  y  accelerometer  Markov  driving  noise  (ft/s5'2) 

1 6.  uBaz  -  z  accelerometer  Markov  driving  noise  (ft/s5'2) 

17.  uBhb  -  baro  altimeter  Markov  driving  noise  (ft/s1/2) 

18.  uB|_j  -  Loran-C  latitude  Markov  driving  noise  (rad/s1/2) 

19.  uBlg  -  Loran-C  longitude  Markov  driving  noise  ( rad/s1'2) 

20.  uSFu  -  Doppler  U  channel  scale  factor  Markov  driving  noise  ( ) 

21.  uBy  -  Doppler  V  channel  boresight  Markov  driving  noise  (rad/s1'2) 

22.  uByy  -  Doppler  W  channel  boresight  Markov  driving  noise  (rad/s1/2) 

23.  uSBn  -  North  sea  bias  Markov  driving  noise  (ft/s3'2) 

24.  uSBe  -  East  sea  bias  Markov  driving  noise  (ft/s3'2) 
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d)  Relationship  Between  Doppler  Velocity  Errors  and  Doppler  Error  States  - 

The  relationship  between  the  actual  errors  in  each  of  the  U,  V  and  W  Doppler  channels  and  the 
Doppler  error  states  being  modelled  in  the  Kalman  filter  can  be  expressed  very  simply  when  only  first- 
order  errors  are  considered.  Define  5U,  5V  and  SW  to  be  the  errors  in  the  U,  V  and  W  Doppler  channels 
respectively.  Let  SFy,  By  and  By;  be  the  Doppler  U  scale  factor  error,  V  boresight  error  and  W  boresight 
error  respectively.  Then  the  relationships  between  these  two  forms  of  the  Doppler  errors  are  as  follows: 

5U  =  SFu U  ;  SV  =  BVU  ;  5W  =  Bw U  (B4) 


where  SFy  is  dimensionless  and  By,  Bw  are  in  the  dimensions  of  radians. 


e)  The  24  Scalar  Differential  Equations  - 

For  each  of  the  24  error  states,  the  explicit  scalar  differential  equation  is  written  out  with  the  non¬ 
zero  terms  ordered  in  such  a  way  that  it  will  be  easy  to  identify  the  corresponding  elements  in  F  and  G. 

8L  =  -  ( vN/Ro2)  5h  +  (1/Ro)&n 

5A.  =  -[VEtanL/(RoCOsL)]5L  +  [vg/(R02oosL)]8h  -  [  1  /  (Rq  cos  L)  ]  Svg 
Sh  =  -  K,  Sh  +  Svz  +  K!  Bhb  +  Kt  uf\, 

SvN  =  [-2|G|veCOsL-ve2/(RoOOs2L)]8L  +  [(vE2tanL  +  vNv2)/R02]5h-(vz/R0)&/N 

-[2|fl|sinL  +  2vEtanL/R0]SvE-(VN/R0)5vz  +  azeE-aEez 
+  Ci  i  Bax  +  C(  2  Bay  +  Ci_3  Baz  +  C)  i  uax  +  C12  nay  +  0^3  uaz  +  ug^ 

Sve  =  [2|Q((vNcosL  +  vzsinL)  +  VEVN /(RoCOs2L)  J8L  +  [  (vE  vz  -  vE  vN  tan  L)  /  Rq2  ]  Sh 
+  (2|fl(sinL+VEtanL/R0)SvN  +  ( VNtanL/Ro -vz/Ro  )&/£ 

-  ( 2  |fl|  cos  L + vE  /  Rq  )  -  az  en  +  3n  ez  +  C2,i  Bax  +  ^2,2  Bay  +  O2  3  Baz 

+  C^i  uax  +  C2,2  u3y  +  C23  uaz  +  ugE 

Svz  =  -  2  |fl|  vE  Sin  L  8L  +  [2  g  /  Rq  -  (ve2  +  Vn2)  /  Rq2  -  K2]  8h  +  2(vN/Ro)&/N 

+  ( 2 |fl| cos L  +  2 vE/ Ro ) SvE  +  a  E  en  -  a  n  ee  -  Sa  +  C31  Bax  +  C3 2 Bay  +  C33  Baz 
+  K2  Bf\,  +  C31  uax  +  C3]2  uay  +  C3  3  uaz  +  ugz  +  K2  uhb 

EN  =  -  |fl|  sin  L  SL  -  vE  /  Ro2  5h  +  (1  /  Rq)  8ve  -  ( |£J|  sin  L  +  vE  tan  L  /  Rq  )  ee  -  (vn/Ro)ez 
+  Cti  Bo*  +  C12  Bwy  +  C1j3  Boz  +  Cij  uwx  +  Ci  ^  uojy  +  Ci^  uwz 

EE  =  (vn  /  Ro2)  Sh  -  (1  /  Rq)  SvN  +  ( |fl|  sin  L  +  ve  tan  L  /  Ro)  cn  -  (|fl|  00s  L  +  vE  /  Ro)  ez 

+  Cj>1  BWx  +  C2,2  Btt)y  +  C2t3  Bo^  +  C2J  UCOx  +  C2,2  U(0y  +  C2,3  U(0Z 

ez  =  [  |Q)  cos  L  +  vE  /  (Rq  cos2L)  ]  SL  -  ve  (tan  L/Rq2)  Sh  +  (tan  L  /  Rq)  SvE 
+  (vn  /  Ro)  en  +  ( |G|  cos  L  +  vE  /  Ro  )  ee  +  03,1  Bo)x  +  C3  2  Ba)y  +  C3  3  Bwz 

+  Cjl  U(l)x  +  C3  2  U(0y  +  C3  3  B(l)Z 

8a  =  K36h  -  K3Bhb  -  K3uT^ 

Bto*  =  -  (1  /tojx)  Ba>„  +  uBox  :  Bcoy  =  -(1 /Toy)  Bcuy  +  uBcay;  Bcoi  =  -  (1 /to,z)  Bo^  +  uBo^ 
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Bax  =  -(1/Xa()Bax  +  uBa„;  Bay  =  -(1  /t^)Bay  +  uBay ;  Ba*  =  -(1 /x®)  Baz  +  uBaz 
Bf\,  =  - (1  /%,)  Bhb  +  uBhb ;  Blat  =  - (1  /tlt)  Beat  +  uBlt;  Blng  =  -  (1  /tlg)Blng  +  uBlg 
SF(j  =  -(1  /Tsf)  SFu  +  uSF(j;  Bv  =  -  (1 /xgv)  By  +  uBy  ;  Bw  =  -  (1  /  xgw)  Bw  +  uBw 
SBN  =  -  (1  / xs)  S8n  +  uSBn  ;  SBg  =  -  (1  / xs)  SBg  +  uSBe 

f)  Identification  of  Non-Zero  Fu’s  ■ 

Based  on  the  explicit  error  state  scalar  differential  equations  shown  above,  the  non-zero  elements 
of  F  can  be  identified  as  follows: 

F  1,3  =  -Vn/Ro2;  Fi  4  =1/Rq 

Fai  =  -VEtanL/(RoCOSL) ;  F^  *  ve/(Ro2cosL)  ;  F2,5  =  -1/(RoCOsL) 
f3,3  =  -  Ki ;  F3  6  =  1  ;  F3 17  =  Ki 
F41  =  -  2  |fi|  vE  00s  L  -  vE2  /  (Ro  cos2  L) 

F43  =  (vE2tan  L  +  vN  vj  /  R<,2  ;  F44  =  -vz/Ro 

F45  =  - 2 |Q| sin L  -  2 vEtan L/ Rq ;  F46  *  -vn/Ro 

F4.8  =  az;  F49  =  -aE;  F414  =  Ct.i  ;  F415  =  1  F416  =  C1>3 

F51  =  2  |fl|  (vN  COS  L  +  vz  sin  L)  +  Ve  vn  /  (Ro  cos2t) 

F53  =  (vevz-ve  vNtanL)/Ro2;  F54  =  2  |flj  sin  L  +  VEtanL/Ro 

F55  =  vNtanL/Ro  -  vz/Ro;  F56  =  -  2  |flj  cos  L  -  vE/Ro ;  F57  =  -az;  F59  -  a^ 

Fs,14  =  C2,i  ;  F515  =  C2,2 :  F5  t6  =  C23 

F6,i  =  -2|fl|vESinL;  F6,3  =  2g/Ro  -  (ve2  +  '^2)/r02  -  k2 

F6,4  =  2 Vn/Ro;  ^6,5  *  2  |flj  COS  L  +  2  VE  /  Ro;  F$,7  =  aE ;  Fgg  =  -aNl  F^io  =  '1 

F6,i4  =  C31  ;  F515  =  C32 ;  F6,i6  =  C33  ;  F6,17  =  K2 

F71  =  -  (0|  sai  L ;  F73  =  -ve/Ro2;  F75  =  1/Rq;  F73  =  -  |fl|  sin  L  -  V£tanL/ Ro 
F79  =  -vn/Ro;  F711  =  C11  ;  F712  =  C12;  F713  =  C13 

F8,3  =  vn/Rq2  ;  Fg_4  =  -1  /  Rq  ;  F87  =  |fl|  sin  L  +  ve  tan  L  /  Ro 
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*  -|Q|cosL  -  vg/Ro;  Fa,n  -  C2,i  ;  Fa,i2  =  C 2,2:  F8,i3  =  02,3 
F9>1  «  |Q|oosL  +  ve/(RoCOs2L);  F9i3  =  -vEtanL/Ro2;  F95  =  tanL/Ro 
F9,7  »  vn/Ro;  F9a  -  |fl(oosL  +  ve/Ro;  F9ih  =  C31 ;  F912  =  03,2;  F913  =  Qjs 
F10,3  =  K3;  F10,t7  *  -K3 

F11,11  =  -1/^xi  F12,12  *  ay;  F13,13  *  -1 

F  14,14  =  -1  /x«  ;  Fi5,i5  =  F^.16  *  -1/Taz 

F17.I7  =  -1  /Thb  ;  F18,18  *  -1  /XLT  ;  F19,19  =  *1/Xuq;  F20,20  =  *1  /  Xg= 

F 21 ,21  =  -1  /Tbv  ;  F22.22  =  -1  /tew;  F 23,23  =  -1  /ts:  F 24,24  =  -1  /Xs 

g)  Identification  of  Non-Zero  Gu’s  - 

The  non-zero  elements  of  G  can  be  identified  as 

G3.10  *  *<1 

G4.4  =  Cii  ;  G45  =  Ci_2 ;  G46  *  C13 ;  G47  *  1 

G5.4  =  C2 1 ;  G5S  =  C22 :  Cs,6  =  02,3 :  G$a  =  i 

G6.4  =  C31 ;  Ge,5  =  Q3.2 ;  Ga6  =  C3i3 ;  G^g  =  1 ;  G6,io  =  K2 

G7.1  =  C1f1  ;  G7S  =  C12  ;  G73  =  C13 

Ga,i  =  0>.i ;  Gas  -  O2.2 :  G*3  =  C2,3 

Gg.i  =  C3,1  ;  G9  2  =  0},2 ;  Gg,3  =  C3i3 

Gio.io  =  -K3;  Gi  1,1 1  =  1;  Gi2,12  =  1;  Gl3,13  *  1  ;  Gi4,i4  =  1 

Gi5,15  =  1  :  Gt6i6  =  1  '•  G1777  =  1  ;  Gia,l8  =  1  ;  Gi9ii9  ■  1 

G20S0  *  1 ;  G21S1  =  1 ;  G2222  *  1 ;  G2323  =  1 ;  G24.24  =  1 
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APPENDIX  C 

Kalman  Filter  Measurement  Processing 


There  are  two  types  of  measurements  to  be  processed  in  the  Kalman  filter:  i)  discrete 
measurements  based  on  position  differences  between  IRS  position  components  and  Loran-C  position 
components  available  every  second;  and  ii)  discrete  measurements  based  on  velocity  differences 
between  IRS  velocity  components  and  Doppler  radar  velocity  components  available  every  half  second.  In 
this  Appendix  the  measurement  error  models  and  observation  equations  are  derived  for  both 
measurement  types. 


a)  Measurements  Based  on  Loran-C  • 

At  any  discrete  time  tk  assume  that  two  simultaneous  measurements  are  processed  by  the  Kalman 
filter  based  on  airborne  Loran-C  position  data,  i.e. 


Zi(k)  =  2l(tk)  =  LATrsOk)  -  LATL0R(tk) 

Z2(k)  =  z2(tk)  *  LNG|Rs(tk)  -  LNGuoR(tk)  (Cl) 


where  LAT  irs,  LNGirs  are  the  IBS  estimates  of  position  and  LAT  lor,  LNGLOr  are  the  Loran-C  estimates 
of  position,  for  the  same  point  in  time  tk.  Let  LAT y^tk)  and  LNGjrOk)  be  the  true  geographical  position  of 
the  aircraft  at  tk-  Define  5LAT|RS,  5LNG|rs  to  be  the  errors  associated  with  the  IRS  calculation  of  position 
and  define  8LATlor.  SLNGlor  to  be  the  errors  associated  with  the  Loran-C  measurement  of  position. 
Then  z^k)  and  z2(k)  can  be  expressed  as 


2l(k)  =.  [LATTr(tk)  +  8LAT,RS(tk)]  -  [LATTf(tk)  +  5UVTL0R(tk)] 
=  aATreflk)  -  5LATLOR(tk) 


Z2(k)  -  (LNGrr(tk)  +  6LNG,Rs(tk)l  -  (LNGrr(tk)  +  8LNGLOR(tk)] 

-  SLNC^tk)  -  8LNGlor(W  (C2) 


A  suitable  error  model  for  the  Loran-C  position  components  would  consist  of  a  bias-like,  first-order 
Markov  process  together  with  additive  ZMWG  random  noise.  Thus,  we  can  write 


5LATLoR<tk)  *  Buvr(tk)  +  vyLT(tk) ;  vllt  ZMWG 
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SLNGLCtt(t|0  =  Bi_NG(tk)  +  Vu_G(tk)  I  Vq.G  ZMWG  (C3) 

where  the  discrete-time  form  of  the  equations  describing  the  Markov  error  processes,  B|_AT(k)  and 
BlngM,  would  be 


BLAT(k+1)  =  exp  (-  AT/tlt)  Bi_AT<k)  +  uBu(k) 

B|_NG(k+1)  =  exp  (-  AT/tlg)  BlnqM  +  uBLo(k)  (C4) 


with  uBlt  and  uBlg  being  the  ZMWG  Markov  driving  noises  associated  with  these  first-order  Markov 
processes.  Note  that  Eqns.  C4  are  simply  the  discrete-time  equivalents  of  the  corresponding  continuous¬ 
time  Markov  equations  for  Loran-C  position  error  states  that  have  been  given  in  Appendix  B. 

The  general  form  of  the  scalar  observation  equation  for  the  ith  measurement  at  tk  is 


zj(k)  =  Hj(k)  8x(k)  +  Vj(k) ;  Vj  ZMWG  (C5) 

and  we  wish  to  identify  the  non-zero  elements  of  Hj  and  H2,  which  will  correspond  to  the  Loran-based 
measurements  that  are  defined  in  Eqn.  Cl.  The  complete  error  state  vector  5x  has  already  been  defined 
in  Appendix  B.  If  we  look  at  Eqns.  C3  in  terms  of  elements  of  5x  we  can  identify 


5LATirs  =  5L ;  5LNGIRS  =  & 

SLATlor  =  Blat  +  VLLT ;  SLNGloR  =  B|_ng  +  Vu_G  (C6) 


Then  Eqns.  C2  can  be  written  in  the  alternative  form 


Zi(k)  =  8L(k)  -  BuT(k)  -  vu.tfk) 

za(k)  =  S^k)  -  BwG(k)  -  vujQ(k)  (C7) 


and  we  can  identify  vi(k)  =  -vu_T(k),  V2(k)  =  -  vu_G(k);  with 

#18 

H,(k)  =  [1  000 . 00-1  000000] 

H2(k)  =[01  0  0 . 0  00-1  00000]  (C8) 

#19 

The  only  non-zero  elements  of  Ht  and  H2  are  seen  to  be  hi  t1  =1,hi18  =  -1 ,  h2,i  =  landtag  =-1. 
The  ZMWG  measurement  noise  components,  vi  and  V2,  are  identified  as  simply  the  negatives  of  the 
corresponding  Loran-C  ZMWG  random  noise  components. 
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b)  Measurements  Based  on  Doppler  Radar  - 

Assume  that  the  fundamental  data  available  from  the  Doppler  radar  system  are  in  the  form  of  a 
three-component  body  axis  velocity  vector  vbo,  with  the  components  denoted  as  Ud,  Vq  and  Wq  (the 
positive  sense  being  forward,  to  the  right  and  down  respectively).  It  is  then  appropriate  to  process  the 
velocity  differences  between  the  IRS  and  the  Doppler  radar  in  the  body  axis  frame.  To  do  this  correctly, 
IRS  velocity  components  have  to  be  transformed  into  equivalent  body  axis  components  and  the  Doppler 
velocity  data  must  be  corrected  for  lever  arm  effects.  The  lever  arm  effects  are  due  to  the  fact  that  the 
Doppler  radar  antenna  is  in  a  different  location  than  that  of  the  IRS.  In  order  to  express  the  observation 
equations  properly,  we  must  also  develop  error  equations  for  each  of  the  IRS  and  Doppler  velocity 
components  in  the  body  axis  frame  of  reference. 

Transformation  of  IRS  Velocity  Vector  into  Body  Axes: 

Let  v<3|  represent  the  three-component  IRS  velocity  vector  in  the  geographic  reference  frame, 
where  the  individual  components  are  vn  ,  vg  and  vz .  Denote  by  C^q  the  transformation  matrix  that  converts 
vGi  into  an  equivalent  velocity  vector  in  body  axis  coordinates  vb|,  having  components  U|,  V|,  W|.  In 
vector/matrix  terminology,  the  transformation  equation  is  simply 


Vb|  =  C^G  V°| 


(C9) 


This  particular  transformation  matrix  can  be  re-computed  continuously  using  Euler  angle  (i.e. 
attitude)  information  from  the  LTN-90-100  IRS.  Note  that  the  transformation  CbG  is  the  inverse 
transformation  of  the  DCM  that  has  already  been  defined  in  Appendix  A  and  denoted  as  CGb-  For 
orthogonal  reference  frames,  such  as  the  ones  we  are  dealing  with,  it  can  be  shown  that 


CbG  -  (CGb)"1  =  (CGb)T 


(CIO) 


Taking  the  transpose  of  CGb,  as  defined  in  Eqns.  All,  results  in  the  following  elements  for  Cbc: 


C  1.1 ' 

= 

cos 0  cosy 

C,,2- 

= 

oos0  siny 

Ci,3' 

= 

sin  6 

C2.1' 

= 

sin  9  sin  $  cos  y  -  cos  $  sin  y 

C  2,2" 

= 

sin  9  sin  $  sin  y  +  cos  <)>  cos  y 

C2.3' 

= 

-  cos  9  sin  <j> 

C3.1 

= 

sin 0  cos 4> cosy  +  sin $  siny 

C  3,2* 

= 

sin  0  cos  <(>  sin  y  -  sin  $  cos  y 

C3,3' 

= 

-  cos  9  COS  «J) 

(C11) 
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Based  on  the  foregoing  definition  of  the  elements  of  CbG>  the  IRS  velocity  components  can  be  expressed 
in  body  axes  as 


Ui 

Vi 

W| 


Cl.l'VN  +  Ci,2'Ve  +  Ct,3'vz 
C2/VN  +  C2^'VE  +  C23'vz 
C3  l'VN  +  C3  2-  VE  +  C3  3'  vz 


(Cl  2) 


Correction  of  Doppler  Velocities  for  Lever  Arm  Effects: 

The  strapdown  IRS  and  the  Doppler  radar  antenna  are  not  co-located  on  the  Twin  Otter  aircraft. 
Due  to  the  length  of  the  lever  arm  between  the  two  locations,  there  will  be  significant  relative  motion 
between  the  two  navaids  caused  by  typical  aircraft  attitude  changes.  It  is  thus  important  to  correct  the 
Doppler  velocity  components  for  this  lever  arm  effect  so  that  the  velocity  difference  measurements, 
referenced  to  the  location  of  the  IRS,  are  not  in  error.  Let  foi  represent  a  three-component  I  ever  arm 
position  vector,  in  body  axis  coordinates,  going  from  the  Doppler  radar  antenna  1c  the  LTN-90-100  IRS. 
The  individual  components  of  Jbi  will  be  defined  as  4,  4,  4  and  they  have  the  values  4  =  13.97  ft,  4  =  0.91 
ft,  lz  =  - 1 .34  ft.  Let  oohg  be  the  angular  body  rate  vector  describing  the  angular  rotation  rate  of  the  IRS 
(and,  hence,  the  aircraft)  as  expressed  in  body  frame  coordinates.  The  components  of  c^b  (i.e.  (dbx,  “By. 
cobz)  are  measured  by  the  IRS,  and  would  be  directly  available  in  digital  form  from  the  IRS  dataport.  Denote 
the  uncorrected  Doppler  velocity  vector  by  vbD  (with  components  Ud.  Vq,  WD)  and  the  corresponding 
corrected  velocity  vector  by  vbDC  (with  components  Udc.  vdc.  wdc)-  The  vector  equation  that  describes 
the  Doppler  velocity  correction  for  leve.  arm  effects  is  then  simply 


VbDC  =  vho  +  fi*>B  X  In 


(Cl  3) 


The  standard  vector  cross  product  formula  results  in 


©bB  x  k i 


r  0  -0& 

“By  1 

m 

“Bz  0 

-“Be 

v 

L'Ofcy  OBx 

0  J 

UJ 

+  (Dq, 

«Sz  4 

-  (0Q< 

4 

L  4 

+  «Bx 

V-l 

(Cl  4) 


The  corrected  Doppler  velocity  components  can  then  be  identified  as 
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Mdc  =  Ud  -  £ifc  V  +  <% 

Vqc  =  Vq  +  Ofe  4  -  COBc  4 

V\fcc  =  Wo  -  Gfy  4  +  WBx  V  (Cl  5) 

Once  the  IRS  velocity  components  have  been  transformed  into  the  body-axis  coordinate  frame 
(Eqns.  C12),  and  the  Doppler  velocity  components  have  been  corrected  for  the  lever  arm  effects  (Eqns. 
Cl 5),  the  three  Doppler-based  measurement  components  can  then  be  computed  as 

Z3(k)  =  U|(k)  -  UDC(k) 

Z4(k)  =  V,(k)  -  VDC(k) 

zs(k)  =  W,(k)  -  WDC(k)  (Cl  6) 

at  any  discrete  time  tR. 


Error  Modelling  for  Doppler-Based  Measurements: 

Let  6vbm  denote  the  three-component  vector  of  differences  between  IRS  and  Doppler  velocities 
as  computed  in  the  body-axis  frame  (the  elements  of  8vbm  are  exactly  the  components  expressed  in 
Eqns.  C16).  A  ■ 'ector/matrix  equation  for  5vbm  would  be 

8vbm  =  vb,  -  vboc 

=  CP G,  VG,  -  vboc  (Cl 7) 


In  order  to  identify  a  suitable  Kalman  filter  observation  equation,  we  would  like  to  derive  an  expression  for 
8vbm  in  terms  of  IRS  and  Doppler  error  quantities.  The  errors  in  CbG,  vQ(  and  vbDC  can  be  expressed  as 

8vG,  =  -  vG  +  (eGCx)vG 

SvbpQ  =  vbpc  -  vb 

SCbG  =  CPo  -  O>0  =  (eGCx)  (Cl 8) 

where 

vG  =  true  velocity  vector  in  geographic  coordinates 
vb  =  true  velocity  vector  in  body  axis  coordinates 
CbQ  =  true  transformation  from  geographic  to  body  axes 
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CbGi  =  transformation  from  geographic  to  body  axes  as  computed  from  error  corrupted  IRS 
attitude  data 

(e00  x)  =  skew  symmetric  matrix  expressing  misalignment  of  IRS  platform  with  respect  to 
computer  frame  (see  Refs.  13, 14) 

Substituting  the  expressions  of  Eqns.  Cl 8  into  Eqn.  Cl 7  results  in 


8vbm  =  0>q  Sv^]  -  CbG  (eGC  x)  (eGC  x)  vG  -  Sv^c  (Cl  9) 

which,  to  first  order,  becomes 

5vbm  =  0>G|  5vG|  -  Svbpc  (C20) 


We  can  identify  the  components  of  8vqi  as 


Svqi  =  [Svn  Sve  Svz]T 


(C21) 


where  8vN,  8ve  and  8vz  are  simply  the  IRS  velocity  error  states  of  the  Kalman  filter  design.  The 
components  of  CbGi  can  be  computed  from  Eqns.  Cl  1,  where  the  attitude  parameters  are  supplied  by  the 
IRS.  It  then  only  remains  to  define  a  suitable  Doppler  error  model  for  8vbD  c  so  that  Svbm  can  be  completely 
specified  in  terms  of  error  states. 

The  main  contributions  to  the  Doppler  velocity  error  vector  8vbDC  will  consist  of  i)  scale  factor 
error,  ii)  boresight  error  (i.e.  Doppler  antenna  misalignment),  iii)  high  frequency  fluctuation  noise  and  iv) 
ocean  current  effects  (when  applicable).  For  a  fixed-wing  aircraft,  in  which  it  can  be  assumed  that  the 
aircraft  velocity  vector  is  approximately  aligned  with  the  U  (longitudinal)  body  axis,  an  adequate  error  model 
for  Svboc  is  the  following: 


r  8u oc  i  r  sfu  i  r  sbn  i  r  nu  i 

SvbQQ  =  5Vqc  =  By  l>r  +  CbQ|  SB  £  +  n  v 

L  5Wqc  J  L  Byy  J  L  0  J  L  nw  J  (C22) 


where 

SFu  =  U  channel  Doppler  scale  factor  error 

By,  Bw  =  V,  W  components  of  Doppler  boresight  error 

Uj  =  true  U  body  velocity  component 

SBn,  SBe  =  north,  east  components  of  sea  bias  (i.e.  the  negative  of  the  actual  ocean  current 
components) 
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nu.  nv ,  nw  =  U,  V,  W  components  of  ZMWG  fluctuation  noise 

The  error  model  of  Eqn.  C22  includes  five  error  states  to  describe  Doppler-based  error 
propagation,  and  these  error  states  are  defined  in  terms  of  the  following  first-order  Markov  processes: 

SFu  =  ( -1  /tsf  )  SFu  +  uSFy  ;  uSFy  ZMWG 
By  =  ( -1  /xgy  )  By  +  uBy  ;  uBy  ZMWG 
Byv  =  ( -1  /lav  )  Bw  +  uByy ;  uBw  ZMWG 
SBjg  =  ( -1  /xs  )  SBfg  +  uSBfg ;  uSBm  ZMWG 

SBe  =  (-1/ts)SBE  +  uSBe  ;  uSBE  ZMWG  (C23) 

These  equations  have  already  appeared  as  part  of  the  set  of  24  scalar  differential  equations  that 
describe  the  complete  set  of  Kalman  filter  error  states  in  Appendix  B.  Substituting  Eqns.  Cll,  C21  and 
C22  into  C20  yields  expressions  for  scalar  measurements  z3l  z4  and  Z5  in  the  general  observation 

equation  form  of  Eqn.  C5.  We  identify  the  following  non-zero  components  of  H3>  H4  and  H5: 

h3  4  =  cos  0  cos  y ;  h3  5  =  cos  0  sin  y ;  h36  =  sin  0 

ha  20  =  -UT;  h323  =  -  cos  0  cosy;  h324  =  -  cos  0  sin  y 

h4  4  =  sin  0  sin  0  cos  y  -  cos  0  sin  y 
h4  5  =  sin  0  sin  $  sin  y  +  cos  0  cos  y;  h4  6  =  -  cos  0  sin  <j> 

h4,2i  =  -  Ut  :  h4  23  =  00s  <}>  sin  y  -  sin  0  sin  0  cos  y 

h4  24  =  -  cos  <}>  cos  y  -  sin  0  sin  0  sin  y 
h5i4  =  sin  0  cos  0  cos  y  +  sin  0  sin  y 
h5  5  =  sin  0  cos  0  sin  y  -  sin  0  cos  y 

h56  =  -cos0cos0;  h522  =  -Ut  (C24) 

The  discrete  measurement  noise  processes  associated  with  z3,  z4  and  Z5  are  identified  as  v3  = 

-  nUt  v4  =  -  ny  and  v5  =  -  nw  respectively.  Thus,  each  measurement  noise  process  has  a  noise  variance 
equal  to  that  of  the  corresponding  Doppler  channel  fluctuation  noise.  Note  that,  for  a  few  of  the  elements 
of  H,  it  is  required  to  know  Ur,  the  true  U  body  velocity  component.  A  reasonable  approximation  to  Ut  can 
be  computed  from 

Ur  =  Uqc  /  ( 1  +  SFy )  (C25) 

where  SFu  would  be  the  latest  value  of  Doppler  scale  factor  error  as  estimated  by  the  Kalman  filter. 
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APPENDIX  D 

Measurement  Averaging  (Prefiltering) 


Recall  that  the  update  rate  for  the  Kalman  filter  has  been  set  at  0.1  Hz,  whereas  the  two 
measurement  types  have  data  available  at  much  higher  data  rates.  In  the  case  of  the  Doppler  velocity 
sensor  the  rate  is  2  Hz,  and  for  the  ARNAV  Loran-C  receiver  the  data  rate  is  1  Hz.  In  order  to  use  all 
measurement  samples  during  the  10  second  Kalman  filter  update  interval,  a  technique  known  as 
measurement  averaging,  or  prefiltering,  is  employed  (see  Ref.  8).  The  use  of  measurement  averaging 
affects  the  computation  of  z,  H,  v  and  R,  as  well  as  the  order  in  which  the  five  Kalman  filter  recursion 
equations  are  executed. 

In  general,  let  AT  =  tk+i  -  tk  be  the  fundamental  Kalman  filter  update  interval,  and  assume  that  a 
set  of  vector  measurements  occurs  at  N  equally  spaced  time  points  within  AT,  i.e. 


z(tk  +  8ti)  ;  i  -  1,2 . N  (D1) 

where  tk  +  Stj  is  the  sample  time  of  the  ith  measurement  within  the  interval  AT.  Recall  the  basic 
vector/matrix  form  of  the  observation  equation: 

z(t)  =  H(t)  8x(t)  +  v(t)  (D2) 


For  the  vector  measurement  at  time  tk  +  8tj  we  can  write 


z(tk  +  5ti)  =  H(tk+8ti)  8x(tk+5ti)  +  v(tk  +  Sti) 


(D3) 


But,  in  general,  we  can  express 


8x(tk  +  8ti)  =  <X>(tk,tk  +  8tj)  8x(tk) 


(D4) 


where  0(tk,  tk  +  8tj)  is  an  intermediate  va^e  of  the  transition  matrix  whose  elements  can  be  expressed  as 
follows: 


tk+Stj 

=  1  Fij(t)dt  ;  i,  j  =  1 . . 24  ;  i*j 

tk 


‘t’i.jOk.  tk  +  Stj) 
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tk+Stj 

«*»i.i(tk.  tk+etj)  =1  +  1  Fj  j(t)  dt  ;  i  =  1, . 10 

tk 

<&i,i(tk.  tk  +  Stj)  =  exp(-8tj/xj )  ;  i  =  11 . 24  (D5) 


Note  that  all  of  the  integrals  in  Eqns.  D5  would  be  computed  numerically,  using  trapezoidal 
integration,  so  that  the  required  intermediate  values  would  automatically  be  available  from  the  trapezoidal 
integration  routines.  Substituting  Eqn.  D4  into  Eqn.  D3  produces 


z(tk  +  Stj)  =  H(tk  +  8tj)  ^tk.tk  +  Stj)  Sx(tk)  +  vdk  +  Stj)  (D6) 


Summing  both  sides  of  Eqn.  D6  over  all  values  of  i  =  1 . N  yields 


N  N  N 

I  z(tk  +  8tj)  =  I  [Hdk+Stj)^,  tk  +  Stj)  8x(tk)]  +  Iv(tk  +  8tj) 
i=1  i=1  i=1 


N  N 

=  [ZH(tk  +  8^)<l»{tk,tk  +  8tj)]  8x(tk)  +  X  v(tk  +  Stj)  (D7) 

i=1  i=1 


Divide  both  sides  of  Eqn.  D7  by  N  to  get 


N  N 

(1/N)  I  z(tk  +  Stj)  =  [  (1/N)  I  H(tk  +  8tj)  4>(tk,  t  k  +  8tj)]  Sx(W 
i=1  i=1 

N 

+  (1/N)  Ivfo+Stj)  '  (D8) 

i=1 

Now  define  the  following  set  of  averaged  quantities,  all  referenced  to  discrete  time  tk: 


N 

z(k)  =  (1/N)  I  z(tk  +  Stj):  averaged  measurement 
i=1 

N 

H(k)  =  (1/N)  I  [H(tk  +  Stj)fc(tk,tk  +  8tj)]:  averaged  measurement  matrix 
i=1 
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N 

v(k)  =  (1/N)  X  v(tk  +  8tj):  averaged  measurement  noise 
i=1 

N 

R(k)  =  (1/N2)  X  R(tk  +8tj):  averaged  measurement  noise  covariance  matrix 
i=1 

=  (1/N)  R;  if  R  is  a  constant  matrix 


With  the  foregoing  definitions,  Eqn.  D8  can  be  written  as 


z(k)  =  H(k)  Sx(k)  +  v(k) 


and  becomes  the  equivalent  averaged  observation  equation,  as  referenced  to  discrete  time  tk- 


(D9) 


COMPUTER  SYSTEM 
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FIG.  1 :  IAR  TWIN  OTTER  ATMOSPHERIC  RESEARCH  AIRCRAFT  AS 
INSTRUMENTED  FOR  AIRBORNE  WIND  MEASUREMENT 
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FIG.  2:  ERROR  STATE  FEEDFORWARD  KALMAN  FILTER  NAVIGATOR 
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FIG.  3:  IRS  /  DOPPLER  /  LORAN-C  KALMAN  FILTER  CONFIGURATION 


NORTH  VELOCITY  (M/S) 
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FIG.  4:  VELOCITY  COMPONENTS  OF  SIMULATED  TRAJECTORY 
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FIG.  5:  VELOCITY  RATE  COMPONENTS  OF  SIMULATED  TRAJECTORY 
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FIG.  6:  ATTITUDE  COMPONENTS  OF  SIMULATED  TRAJECTORY 
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FIG.  7:  ATTITUDE  RATE  COMPONENTS  OF  SIMULATED  TRAJECTORY 
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FIG.  8:  POSITION  CHANGE  COMPONENTS  OF  SIMULATED  TRAJECTORY 


LATITUDE  ERROR  (NM) 


10  20  30  40 

LONGITUDE  ERROR  (NM) 


50 


-i ■ -  •  —  r  . . r  — . — — - - - - - 1 - - - - - - - j- 

'0  20  30  *0  50 

ALTITUDE  ERROR  (FT) 


TIME  (minutes) 

FIG.  9:  STRAPDOWN  NAVIGATOR  POSITION  ERRORS 
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FIG.  10:  STRAPDOWN  NAVIGATOR  VELOCITY  ERRORS 
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FIG.  1 1 :  STRAPDOWN  NAVIGATOR  ATTITUDE  ERRORS 
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FIG.  12:  LORAN-C  LATITUDE  AND  LONGITUDE  MARKOV  PROCESSES 
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FIG.  13:  DOPPLER  VELOCITY  MARKOV  PROCESSES 


0000  -  NO.  OF  SECONDS  OF  INITIAL  DATA  TO  BE  SKIPPED 

03590  -  DURATION  OF  KF  RUN  (IN  SECONDS,  DIVISIBLE  BY  10) 

00000  -  START  TIME  FOR  OVERWATER  FLIGHT  (IN  SECS  RELATIVE) 

00000  -  END  TIME  FOR  OVERWATER  FLIGHT  (TO  START  OF  FILE) 

82.500  -  NOMINAL  OPERATING  LATITUDE  ( DEGS  ) 
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FIG.  14:  LISTING  OF  FILE  SDKFNAV  CTRL 
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FIG.  15:  KALMAN  FILTER  LORAN-C  POSITION  MEASUREMENTS 
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FIG.  16:  KALMAN  FILTER  DOPPLER  VELOCITY  MEASUREMENTS 
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FIG.  17:  ACCURACY  OF  IRS  POSITION  ERROR  STATE  ESTIMATION 
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FIG.  18:  ACCURACY  OF  IRS  VELOCITY  ERROR  STATE  ESTIMATION 
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FIG.  19:  IRS  TILT  AND  VERTICAL  ACCELERATION  ERROR  STATE  ESTIMATES 
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FIG.  20:  IRS  GYRO  DRIFT  ERROR  STATE  ESTIMATES 
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FIG.  21 :  IRS  ACCELEROMETER  AND  BAROMETRIC 
ALTIMETER  ERROR  STATE  ESTIMATES 
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FIG.  22:  ACCURACY  OF  LORAN-C  POSITION  AND  DOPPLER  SEA  BIAS 

ERROR  STATE  ESTIMATION 
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FIG.  25:  INNOVATIONS  FOR  KALMAN  FILTER  VELOCITY  MEASUREMENTS 
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FIG.  26:  HORIZONTAL  TRACK  PLOT  FOR  FLIGHT  #2 
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FIG.  27:  LORAN  -C  POSITION  MEASUREMENTS  FOR  FLIGHT  #1 
LORAN-C  DATA  PROBLEMS 
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FIG.  28:  LORAN-C  POSITION  MEASUREMENTS  FOR  FLIGHT  #1 
PROBLEMS  CORRECTED 
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FIG.  29:  DOPPLER  VELOCITY  MEASUREMENTS  FOR  FLIGHT  #1 
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FIG.  30:  IRS  POSITION  ERRORS  AT  VOT'S  FOR  FLIGHT  #1 


0.0  to.o  20.0  30.0  40.0  50.0  60.0  70.0  80.0  90.0 

TIME!  MINS.) 


FIG.  31 :  LORAN-C  POSITION  ERRORS  AT  VOT’S  FOR  FLIGHT  #1 
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0000  -  NO.  OF  SECONDS  OF  INITIAL  DATA  TO  BE  SKIPPED 

05210  -  DURATION  OF  KF  RUN  (IN  SECONDS,  DIVISIBLE  BY  10) 

00000  -  START  TIME  FOR  OVERNATER  FLIGHT  (IN  SECS  RELATIVE) 

00000  -  END  TIME  FOR  OVERWATER  FLIGHT  (TO  START  OF  FILE) 

45.000  -  NOMINAL  OPERATING  LATITUDE  (DEGS) 

-06.63  000.92  001.48  -  U,V,W  ACCELER-TO-INS  LEVER  ARMS  (FEET) 

*  INS  SYSTEM  ERROR  STATE  STATISTICS  (1-SIGMA): 

I.  0 . 2000E00  2.  O.ZOOOEOO  LAT  8  LONG  (NM) 

3.  3 . 0000E01  4.  0.3300E01  ALT  (FT)  &  N.  VELOC.  (FT/SEC) 

5.  0.3300E01  6.  0.3300E01  EAST  8  VERT.  VELOC'S  (FT/SEC) 

7.  0.3340E-1  8.  0.3340E-1  N.  8  E.  AXIS  ATTITUDE  (DEG) 

9  .  0.8594E-1  10.  0.2500E-1  Z  AXIS  ATT.  (DEG)  8  BAP.O  LOOP  ACC.  ( FT/S*#2 ) 

*  SENSOR  MARKOV  ERROR  STATISTICS  (1-SIGMA  VALUE  8  CORR.  TIME  IN  SECS.): 

II.  2 . 0G00E-2  07500.0  X  AXIS  GYRO  (DEG/HR) 

12.  2.0000E-2  07500.0  Y  AXIS  GYRO  (DEG/HR) 

13.  2 . 0000E-2  07500.0  Z  AXIS  GYRO  (DEG/HR) 

14.  1 . 0000E-2  07500.0  X  AXIS  ACCEL  ( FT/SEC*#2  ) 

15.  2.000CE-3  15000.0  Y  AXIS  ACCEL  ( FT/SEC**2 ) 

16.  2 . OOOOc-3  15000.0  Z  AXIS  ACCEL  ( FT/SEC**2  ) 

17.  2.0000E02  05000.0  ALTIMETER  BIAS  (FEET) 

18.  0 . 2000E-0  03600.0  LATITUDE  LORAN-C  FIX  (NM) 

19.  0. 2000E-0  03600.0  LONGITUDE  LORAN-C  FIX  (NM) 

20.  0 . 0200E00  03600.0  DOPPLER  SCALE  FACTOR 

21.  2 . 3000E00  03600.0  OOPPLER  V  BORESIGHT  (DEG) 

22.  1. 1500E00  03600.0  DOPPLER  V)  BORESIGHT  (DEG) 

23.  4.5QOOEOO  00900.0  NCRTH  SEA  BIAS  (FT/SEC) 

24.  4.5000E00  00900.0  EAST  SEA  BIAS  (FT/SEC) 

*  INS  SYSTEM  PLANT  NOISE  STATISTICS  (1-SIGMA): 


1.  2. 

3000E-6  2. 

2. 3000E-6 

INTEGRATED  X  8  Y  GYRO  (DEG) 

3.  2 

3000E-6  4. 

1 . 0000E-3 

INTEGRATED  Z  GYRO  (OEG)  8  X  ACC.  (FT/SEC) 

5.  1 

0000E-3  6. 

1 . 0000E-3 

INTEGRATED  Y  8  Z  ACCEL.  (FT/SEC) 

7.  1. 

0000E-3  8. 

1 . 0000E-3 

INTEGRATED  N  8  E  GRAVITY  (FT/SEC) 

9.  1 

0000E-3  10. 

3 . 1620E01 

INTEG'D  Z  GRAVITY  (FT/S)  8  ALT  NOISE  (FT-S) 

*  MEASUREMENT  NOISE  STATISTICS  (UNAVERAGED,  1-SIGMA): 

1.  0 

0700E00  2. 

O.O7O0EOO 

LAT  8  LONG  LORAN-C  (NM) 

3.  5 

3400E00  4. 

1 . 0680E01 

5.  5.3400EOO  U,V,H  DOPPLER  (FT/SEC) 

6.  5 

OOOOEOO 

BAROMETRIC  ALTIMETER  (FEET) 

»  VISUAL  ON-TOP  REFERENCE  INFORMATION: 

0  15 

-  POS'N  FIX  FLAG  8  TOTAL  NO.  OF  POINTS 

l.OOOE-l  l.OCOE-1 

-  1  SIGMA  LEVEL  (NM)  FOR  LAT  8  LONG  FIX 

0490 

45.44135 

75.89501 

1  -  TIME,  LAT/LONG  8  FLAG  FOR  EACH  ON-TOP 

0910 

45.47167 

76.23333 

1 

1470 

45.58434 

76.67900 

1 

2010 

45.31082 

76.89049 

1 

2270 

45.19934 

76.90734 

1 

2490 

45.24216 

76.76033 

1 

2840 

45.38618 

76.57565 

1 

3520 

45.18333 

76.12667 

1 

3960 

45.07417 

75.90467 
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45.16051 

75.83766 
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45.16008 

75.83814 

1 

4400 

45.03076 

75.90078 

1 

4540 

45.07919 

75.90170 

1 

4710 

45.15884 

75.83833 
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5200 
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FIG.  32:  SDKFNAV  CTRL  FILE  FOR  FLIGHT  #1 
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FIG.  33:  IRS  POSITION  ERROR  STATE  ESTIMATES  FOR  FLIGHT  #1 


TIME  (minutes) 


FIG.  34:  IRS  VELOCITY  ERROR  STATE  ESTIMATES  FOR  FLIGHT  #1 
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FIG.  35:  IRS  TILT  AND  VERTICAL  ACCELERATION  ERROR 
STATE  ESTIMATES  FOR  FLIGHT  #1 
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FIG.  36:  IRS  GYRO  DRIFT  ERROR  STATE  ESTIMATES  FOR  FLIGHT  #1 
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FIG.  37:  IRS  ACCELEROMETER  AND  BAROMETRIC  ALTIMETER 
ERROR  STATE  ESTIMATES  FOR  FLIGHT  #1 
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FIG.  38:  LORAN-C  POSITION  AND  DOPPLER  SEA  BIAS  ERROR  STATE 

ESTIMATES  FOR  FLIGHT  #1 
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FIG.  39:  DOPPLER  VELOCITY  ERROR  STATE  ESTIMATES  FOR  FLIGHT  #1 
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FIG.  40:  KALMAN-CORRECTED  IRS  POSITION 
ERRORS  AT  VOT'S  FOR  FLIGHT  #1 
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FIG.  40a:  INNOVATIONS  FOR  POSITION  MEASUREMENTS  FROM  FLIGHT  #1 
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FIG.  41 :  LORAN-C  POSITION  MEASUREMENTS  FOR  FLIGHT  #2 
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FLIGHT  #  2  -  LONGITUDE  ERROR  (NM) 
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FIG.  43:  IRS  POSITION  ERRORS  AT  VOT'S  FOR  FLIGHT  #2 
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FIG.  44:  LORAN-C  POSITION  ERRORS  AT  VOT'S  FOR  FLIGHT  #2 
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FIG.  45:  IRS  POSITION  ERROR  STATE  ESTIMATES  FOR  FLIGHT  #2 
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FIG.  46:  IRS  VELOCITY  ERROR  STATE  ESTIMATES  FOR  FLIGHT  #2 
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FIG.  47:  IRS  TILT  AND  VERTICAL  ACCELERATION  ERROR 
STATE  ESTIMATES  FOR  FLIGHT  #2 
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FIG.  48:  IRS  GYRO  DRIFT  ERROR  STATE  ESTIMATES  FOR  FLIGHT  #2 
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FIG.  49:  IRS  ACCELEROMETER  AND  BARO  ALTIMETER  ERROR 
STATE  ESTIMATES  FOR  FLIGHT  #2 
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FIG.  50:  LORAN-C  POSITION  ERROR 
STATE  ESTIMATES  FOR  FLIGHT  #  2 
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FIG.  51 :  DOPPLER  VELOCITY  ERROR  STATE  ESTIMATES  FOR  FLIGHT  #2 
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FIG.  52:  KALMAN-CORRECTED  IRS  POSITION  ERRORS 
AT  VOT'S  FOR  FLIGHT  #2 
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FIG.  53:  COMPARISON  OF  KALMAN  FILTER  CONFIGURATIONS  - 
LATITUDE  ERRORS  FOR  FLIGHT  #1 


FIG.  54:  COMPARISON  OF  KALMAN  FILTER  CONFIGURATIONS  - 
LONGITUDE  ERRORS  FOR  FLIGHT  #1 
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FIG.  55:  COMPARISON  OF  KALMAN  FILTER  CONFIGURATIONS 
LATITUDE  ERRORS  FOR  FLIGHT  #2 
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FIG.  56:  COMPARISON  OF  KALMAN  FILTER  CONFIGURATIONS  - 
LONGITUDE  ERRORS  FOR  FLIGHT  #2 
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FIG.  57:  COMPARISON  OF  KALMAN  FILTER  CONFIGURATIONS  - 
NORTH  VELOCITY  ERRORS  FOR  FLIGHT  #1 
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FIG.  58:  COMPARISON  OF  KALMAN  FILTER  CONFIGURATIONS 
EAST  VELOCITY  ERRORS  FOR  FLIGHT  #1 
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FIG.  60:  COMPARISON  OF  KALMAN  FILTER  CONFIGURATIONS  - 
EAST  VELOCITY  ERRORS  FOR  FLIGHT  #2 
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